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INTRODUCTION 


The  American  Cancer  Society  estimated  that  there  would  be  approximately  230,110  new  cases 
diagnosed  and  approximately  29,900  prostate  cancer  related  deaths  in  2004  [1]  Prostate  cancer 
screening  today  generally  uses  the  Prostate  Specific  Antigen  (PSA)  blood  testing,  free  PSA  testing 
and  Digital  Rectal  Examination  (DRE).  When  using  a  ‘cutoff’  of  PSA  >  4.0  ng/mL  and  an  abnormal 
DRE,  sensitivity,  specificity  and  Positive  Predictive  Value  (PPV)  are  38%,  88%  and  56% 
respectively  [2],  When  either  an  elevated  PSA  or  an  abnormal  DRE  are  used,  (in  isolation  -  not  in 
combination),  sensitivity,  specificity  and  PPV  are  even  lower  [2],  When  the  PSA  is  used  there  exists 
a  significant  gray  area  (4-10  ng/mL)  in  which  cancers  may  be  missed  and  yet  the  number  of 
negative  biopsies  is  large.  Even  though  cancer  detection  sensitivity,  specificity  and  PPV  are 
improved  by  combining  PSA  and  DRE  [2,  3]  the  usefulness  of  DRE  remains  fundamentally  limited 
due  to  its  subjective  nature.  Additionally,  DRE  is  practically  limited  to  the  detection  of  shallow 
(subcapsular)  palpable  abnormalities.  Even  systematic  multi-core  biopsy  fails  to  detect  clinically 
detectable  cancers  in  up  to  34%  of  men  [4].  However,  there  is  evidence  that  as  additional  biopsies 
cores  are  added,  sensitivity  improves  [5],  This  observation  has  resulted  in  an  increase  the  number 
of  cores  taken  during  routine  examination.  Nevertheless,  biopsy-based  detection  sensitivity  remains 
less  than  ideal.  Thus,  there  is  plenty  of  compelling  clinical  interest  in  finding  improved  methods  for 
the  early  diagnosis  of  prostate  cancer  with  improved  sensitivity  and  specificity.  One  recent  example 
of  progress  in  the  field  of  prostate  cancer  detection  involves  an  effort  to  automate  the  DRE 
examination.  Savazyan  recently  described  a  system  for  'mechanical  imaging'  of  the  prostate  [6]. 

This  system  comprises  a  rectal  probe  that  is  instrumented  with  an  array  of  pressure  sensing  strain 
gages  and  a  3D  magnetic  positioner  device.  In  an  in  vitro  trial  [7],  the  new  system  correctly  detected 
and  located  100%  of  the  nodules  under  examination.  This  compares  with  detection  rates  of  83% 
and  67%  for  an  experienced  urologist  and  a  student  respectively.  Thus,  a  significant  improvement 
over  the  conventional  DRE  examination  has  been  demonstrated  for  the  in  vitro  case.  Another  recent 
development  is  the  observation  that  the  sensitivity  of  an  ultrasound  examination  can  be  improved  by 
the  use  of  a  microbubble  based  contrast  agents  [4],  Frauscher's  approach  [4]  involved  the  use  of 
contrast  agent  enhanced  Color  Doppler  that  improved  the  detection  of  hypervascular  regions 
associated  cancer.  Prostate  cancer  was  detected  by  contrast  agent  assisted  ultrasound  in  23  of  24 
patients  known  to  have  prostate  cancer.  (The  method  used  for  determining  definitively  which 
patients  had  cancer  is  not  entirely  clear  in  the  article.)  In  comparison,  conventional  ultrasound 
detected  cancer  in  17  patients.  The  contrast  agent  assisted  approach  detected  cancer  in  8  patients 
with  a  negative  systematic  biopsy-based  diagnosis.  However,  the  cost  of  the  contrast  agent  used  in 
this  study  was  $65  per  patient.  This  cost  makes  up  approximately  half  of  the  cost  of  a  conventional 
ultrasound  examination  and  therefore  represents  a  considerable  impediment  to  its  widespread 
acceptance.  However,  more  recent  publications  [8,  9]  (including  one  from  Frauscher's  group)  cast 
doubt  on  the  true  extent  of  the  improvement  in  diagnostic  accuracy  obtained  by  using  contrast 
agents.  Specifically,  Halpern  was  unable  to  detect  cancers  in  the  inner  gland  and  achieved  a  cancer 
detection  sensitivity  of  only  42%  [8], 
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BODY 


The  work  conducted  as  part  of  this  Army  funded  program  can  be  considered  as  divided  into  the  following  key 
“Aims”: 

1 .  Research,  design,  development  and  prototype  testing  of  a  new  transrectal  ultrasound  transducer, 
syringe  pump  and  ultrasound  instrumentation  to  facilitate  a  Synthetic  Digital  Rectal  Examination 
(SDRE). 

2.  Research,  development  and  prototype  testing  of  techniques  to  enable  quantitative  (dimensionally 
accurate)  3D  reconstructions  of  the  prostate 

3.  Research,  development  and  test  of  techniques  to  improve  ultrasound  image  quality  and  to  facilitate 
automated  (or  semi-automated)  border  detection  of  lesions 

4.  Small  scale  clinical  test  at  the  University  of  Virginia 

Progress  has  been  made  in  each  of  these  tasks  in  the  second  year  of  the  grant.  The  work  is  on-track. 

(The  fourth  section  of  work,  the  small  scale  clinical  test,  will  occur  in  the  third  year  of  the  grant.) 

Progress  with  respect  to  the  areas  are  related  directly  to  the  committed  Statement  of  Work  that  was  funded: 

Aim  1.  Design,  specify,  and  have  built,  a  high  resolution  transducer  optimized  for  imaging 
elastic  inhomogeneities,  unsurpassed  B-Mode  image  resolution  and  possessing  integrated 
3D  capability. 

A  high  frequency  (8-1 4  MHz)  transducer  array  was  designed  and  specified  as  committed  in  Year  1 . 
Vermon  SA,  Tours,  France  made  the  transducer  and  it  was  delivered  in  Q1  2005.  The  transducer 
has  two  tracking  arrays  each  with  32  elements  and  a  central  imaging  array  with  192  elements.  The 
elements  are  spaced  on  a  0.2  mm  pitch.  This  transducer  is  providing  unsurpassed  imaging 
resolution  in  a  transducer  housing  designed  for  transrectal  ultrasound.  The  transducer  will  provide 
the  very  best  image  data  as  a  solid  foundation  for  the  subsequent  work  elements.  This  transducer  is 
operable  at  up  to  14MHz  whereas  the  previously  available  transducer  was  only  operable  up  to  8 
MHz.  Consequently,  we  are  observing  excellent,  and  significant,  improvements  in  image  resolution. 

Aim  2.  Develop  and  test  a  tissue  elasticity  imaging  system. 

As  committed,  we  have  assembled  the  apparatus  to  enable  the  new  approach  to  transrectal 
ultrasound  based  strain  imaging.  (Most  of  this  work  was  completed  in  Year  1.)  We  have  also 
fabricated  several  custom  prostate  phantoms  using  locally  developed  techniques  [10],  By  making 
the  phantoms  internally,  we  are  able  to  iterate  efficiently  the  design  and  also  to  fabricate 
replacement  phantoms  at  very  low  cost  in  a  timely  manner.  Phantoms  tend  to  deteriorate  over  time 
due  to  dehydration  through  the  membranes.  We  have  tested  the  tissue  elasticity  system  using  both 
an  older  8  MHz  transrectal  transducer  and  the  newer  14  MHz  transducer  connected  to  our  Siemens 
Sequoia  ultrasound  machine.  Our  techniques  can  be  migrated  to  other  ultrasound  systems  if 
resource  and  contractual  issues  are  addressed. 

The  deliverable  for  this  phase  is  complete  as  of  the  end  of  Year  2  work  and  a  validated  (phantom) 
elasticity  imaging  system  using  the  approach  presented  in  the  proposal  has  been  produced.  We 
have  verified  our  ability  to  obtain  very  fine  resolution  3D  reconstructions  and  fine  resolution  elasticity 
images  using  the  newly  completed  tissue  elasticity  system. 
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Aim  3.  Apply  image  processing  algorithms  to  acquired  B-Mode  and  elasticity  images  to 
improve  the  quantification  of  detected  elastic  inhomogeneities 


We  have  developed  3D  surface  rendering  from  2D  slices  by  implementing  a  3D  gradient  vector  flow 
(GVF)  snake  algorithm  [1 1],  The  algorithm  relies  on  the  edges  to  be  well  delineated  and  contrast 
between  the  various  regions  to  be  well  defined  in  each  2D  ultrasound  slice  to  produce  a  surface  that 
resembles  the  actual  scanned  object.  A  preprocessing  despeckling  step  is  needed  to  reduce  the 
variances  in  pixel  values  within  homogeneous  regions.  We  have  evaluated  a  wide  variety  of  well 
known  methods  and  a  novel  stochastically  driven  method  design  specifically  for  3D  surface 
rendering  from  2D  slices  of  the  prostate  and  other  organs.  The  method  we  develop  specifically  for 
the  task  of  this  grant  is  a  stochastically  driven  compression  filter  called  the  squeeze  box  filter  (SBF). 
Our  quantitative  evaluation  using  a  Field  II  [12]  simulated  B  mode  ultrasound  image  with  contrast 
enhancement  performance  determined  by  a  modified  Fisher  discriminant  has  determined  that  the 
newly  developed  SBF  method  outperformed  the  other  methods  and  is  exceptional  in  providing  the 
needed  intraregion  reduction  in  variance  and  inter-region  contrast  enhancement  with  computational 
efficiency. 

Excellent  progress  with  respect  to  Aim  3  has  been  made  in  the  second  year  of  the  grant. 

The  fourth  Aim  from  the  Statement  of  Work  relates  to  a  small  scale  clinical  validation  in 
collaboration  with  partially  funded  University  of  Virginia  collaborator  -  Dr.  Dan  Theodorescu.  This 
item  of  work  will  be  addressed  in  the  final  year  of  the  grant. 
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Summary  Statement  of  work  completed  to  date: 

(Comprises  part  of  work  included  as  Appendix  plus  recent  image  processing  results.) 

Prostate  Phantom 

A  purpose-built  prostate  phantom  was  designed  using  the  method  described  in  Negron  et  al.  [10]. 
A  simulated  lesion  (approximately  0.3  ml)  was  formed  inside  of  an  egg-shaped  tissue  region 
mimicking  a  prostate.  A  cylindrical  cavity  (20  mm  in  diameter)  was  formed  to  mimic  the  anal 
opening  and  to  allow  access  for  the  transrectal  transducer.  A  hypoechoic  gel  component 
surrounds  these  three  components.  (Strictly  speaking,  the  surrounding  gel  should  be  echogenic 
too  but  the  lack  of  echoic  inclusions  is  immaterial  in  our  phantom  studies.)  The  lesion  is  made  of 
17%  (by  weight)  acrylamide  gel;  the  tissue  and  exterior  component  are  made  of  5%  acrylamide 
gel.  Thus,  the  lesion  is  perceptibly  stiffer  (approximately  10  times  stiffer )  than  both  the  egg 
shaped  tissue  region  and  the  exterior  gel  component.  A  similar  concentration  (by  volume)  of 
Sephadex  was  added  to  both  the  lesion  and  the  tissue,  resulting  in  similar  ultrasound  image 
intensity  in  these  two  structures  (The  lesion  was  made  slightly  brighter  than  the  tissue  in  order  to 
assist  navigation  during  scanning.)  A  B-mode  ultrasound  image  of  the  prostate  is  shown  in  Fig.  1 . 
The  lesion  in  the  image  is  almost  isoechoic.  This  phantom  is  similar  to  ones  we  have  made  since 
the  beginning  of  the  project. 


Legend: 

1  -  lesion  mimicking 
inclusion 

2  -  prostate  tissue 
mimicking 
surrounding 

3  -  rectal  cavity 

4  -  exterior  casing. 


Fig.  1 .  Left:  Schematic  of  the  transrectal  design  of  prostate  phantom.  Right:  An 
ultrasound  image  of  the  phantom. 


Transducer 

Vermon  SA,  Tours,  France  fabricated  an  8-14MHz  transrectal  transducer  to  our  specification.  This 
was  delivered  in  Q1  2005.  Most  prostate  transducers  used  today  in  premium  ultrasound  scanners 
use  a  tightly  curved  array  placed  on  the  end  of  the  transrectal  probe.  Thus,  these  probes  have 
limited  aperture  contribution  to  image  resolution  at  any  particular  point  in  the  image  field.  However, 
this  new  transducer  has  a  linear  array  format  and  hence  has  a  longer  available  aperture  that  results 
in  finer  lateral  resolution.  Thus,  we  believe  that  this  transducer’s  imaging  resolution  is  practically 
unmatched  in  the  field  of  prostate  ultrasound.  Image  resolution  is  approximately  0.2  mm  lateral  and 
0.1  mm  in  the  axial  (range)  dimension.  The  array  pitch  is  0.2  mm.  There  are  192  elements  in  the 
imaging  array  and  32  elements  in  each  of  the  two  perpendicular  tracking  arrays  that  provide  the 
transducer  with  “I-Beam”  3D  tracking  [13].  This  form  of  3D  tracking  yields  approximately  4.6% 
accuracy  at  the  two  standard  deviation  level.  (95%  of  measurements  will  be  within  4.6%.)  The  I- 
Beam  transducer  is  also  uniquely  matched  to  the  transrectal  prostate  ultrasound  application  for  the 
following  reasons:  a)  the  transrectal  probe  with  the  tracking  mechanism  near  the  transducer 
minimizes  numerical  ill-conditioning  that  may  arise  if  the  means  of  tracking  is  separated  from  the 
imaging  array,  and  b)  the  I-Beam  transducer  estimates  the  relative  tissue  motion  rather  than 
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absolute  tissue  motion,  which  enables  efficient  and  accurate  measurements  even  if  there  is  patient 
motion  of  the  type  that  can  defeat  a  3D  system  that  uses  a  fixed  origin  for  3D  positioning  (eg. 
Magnetic-based,  optical-based,  or  articulated  arm-based  positioner.) 

The  newly  developed  transducer  required  new  software  drivers  to  be  developed  to  enable  it  to  “run” 
on  our  Siemens  Sequoia  research  scanner.  One  graduate  research  assistant  (Yinbo  Li)  spent  the 
summer  of  2005  at  Siemens  Engineering  in  Mountain  View,  CA,  developing  the  required  software 
with  assistance  from  Siemens  Engineering  Staff  -  primarily  Greg  Holley.  Siemens  assisted  with 
costs  associated  with  this  development.  Siemens  Engineering  also  performed  ultrasound  output 
intensity  measurements  to  verify  that  the  transducer  satisfies  current  FDA  regulated  intensities 
(primarily  that  Mechanical  Index  (Ml)  <=  1 .9).  While  the  transducer  was  having  software 
development,  we  also  took  advantage  of  the  opportunity  to  add  a  contrast  agent  imaging  mode  for 
potential  future  work  in  this  area.  Contrast  Pulse  Sequences  [14]  was  implemented  on  the 
transducer.  It  is  intended  that  this  will  enable  future  work  that  might  be  based  on  measuring 
perfusion  in  prostate  or  locating  the  presence  of  molecular  targeted  ultrasound  contrast.  These 
applications  are  beyond  the  scope  of  the  currently  funded  work  and  will  not  be  pursued  without 
future  funding  and  any  requisite  permission. 

The  transducer,  system  and  phantom  are  assembled  into  a  complete  working  3D  scanning  / 
elastographic  system  by  adding  a  latex  sheath  over  the  transducer  (secured  with  elastic  bands), 
Tygon™  flexible  plastic  tubing  and  syringe  pump  to  controllable  inflate  the  sheath  with  plain  tap 
water.  When  these  components  are  assembled  we  have  the  basic  apparatus  for  the  “Synthetic 
Digital  Rectal  Examination”  described  in  the  proposal.  The  programmable  syringe  pump  is  a 
Harvard  Instruments  PHD  2000,  (Harvard  Apparatus,  Holliston  MA).  This  pump  enables  automatic 
water  inflation  and  can  generate  a  quasistatic  stress  and  produce  as  uniform  tissue  deformation  as 
possible.  A  syringe  volume  of  60  ml  was  chosen  to  provide  sufficient  water  to  compress  and  deform 
the  rectal  wall  thus  providing  an  optimal  tissue  strain.  This  volume  is  also  appropriate  in  that  when 
used  in  a  clinical  setting,  the  ultimate  size  of  the  syringe  makes  the  water  injection  process  safe  in 
that  the  syringe  is  emptied  before  any  patient  injury  could  be  anticipated.  We  have  recently 
discovered  that  other  research  groups  have  also  adopted  a  somewhat  similar  balloon  inflation 
method  but  that  these  earlier  efforts  use  a  manually  operated  syringe  [15,  16]. 


Fig.  2.  The  transrectal  transducer  is  covered  with  a  latex  condom.  Water  was  inflated  by 
the  syringe  during  imaging 


Fig.  3.  I-Beam  transducer  -  possessing  a  main  imaging  array  in  the  center  and  a 
tracking  array  in  each  end 

Five  elevational  slices  each  comprising  of  100  Sequoia  14  MHz  In-phase/Quadrature  (l/Q)  data 
frames  were  acquired  during  the  in-vitro  experiment  using  an  incremental  strain  between  the 
consecutive  frames  of  0.04  %.  Total  applied  strain  between  the  first  and  the  last  frame  was 
approximately  4  %.  This  strain  is  broadly  in  accordance  with  the  degree  of  strain  that  has  been  found 
to  be  optimal  for  strain  imaging  [17].An  example  B-mode  image  obtained  by  the  transrectal  I-Beam 
transducer  is  shown  in  Fig.  4. 


Fig.  4.  A  B-mode  image  acquired  with  the  I-Beam  transrectal  transducer  from  the 
prostate-mimicking  phantom,  showing  the  layout  of  image  planes  formed  by 
‘Tracking’  arrays  and  ‘Imaging’  array.  T  -  tracking  array,  T  -  main  imaging  array. 
Delineated  is  a  lesion-mimicking  inclusion. 


Elasticity  Calculation  and  3D  Reconstruction 

The  acquired  IQ  data  were  filtered  using  a  low-pass  filter  to  reduce  jitter,  electronic  noise,  and  out-of 
band  noise.  Six  pairs  of  frames  with  differential  strain  of  2%  were  tracked  using  a  time-domain 
cross-correlation  technique.  Signal  “companding”  or  stretching  was  employed  to  maximize  the 
cross-correlation  coefficient  between  the  pre-  and  the  post-compression  frames.  Companding 
techniques  improve  contrast  to  noise  ratio  of  the  strain  images.  This  improvement  is  desirable  as 
the  prostate  lesions  are  known  to  be  twice  or  at  most  thrice  stiffer  than  the  prostate.  Sub-sample 
precision  was  obtained  in  the  delay  estimates  by  using  a  quadratic  fit  to  the  cross-correlation 
function.  A  search  window  of  approximately  5  wavelengths  was  used  for  time-delay  estimation. 
Lateral  motion  was  tracked  using  the  technique  described  by  Lubinski  et  al.[18].  Displacement 
estimated  from  these  six  different  renditions  were  averaged  to  eliminate  uncorrelated  random  noise. 
Strain  estimates  were  then  obtained  by  taking  the  local  gradient  of  the  displacement  image  (in  axial 
direction).  The  above  process  was  repeated  by  sweeping  the  transducer  in  the  elevational  direction 
and  six  elevational  slices  were  obtained.  In  the  elevation  direction,  the  inter-slice  distance  was 
estimated  with  a  block  matching  approach  based  upon  the  minimum  sum  of  absolute  differences 
(MSAD)  algorithm  of  the  I-Beam  ‘Tracking’  data.  The  inter-slice  distance  was  found  to  be  slightly 
larger  in  the  deeper  portion  than  in  the  shallower  portion,  indicating  that  the  transducer  was  rotated 
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by  small  angles  between  measurements.  In  the  measurements,  the  slices  were  sampled  with  an 
interval  no  greater  than  2  mm,  and  the  rotation  angle  increments  were  less  than  2  degrees  in  order 
to  ensure  that  the  elevational  motion  could  be  accurately  tracked.  Once  the  elevational  motion  of 
these  blocks  was  calculated,  the  acquired  image  slices  were  interpolated  on  to  a  regular  3D  grid  in 
Matlab  enabling  3D  volumes  to  be  rendered 

Figure  5  illustrates  the  lesion  detection  process.  The  elasticity  reconstruction  algorithm  was  tested 
on  the  phantom  data.  The  elasticity  reconstruction  process  for  the  detection  of  prostate  cancer  is 
complicated  than  the  elasticity  reconstruction  process  of  the  breast  cancers.  Unlike  the  boundary 
conditions  for  breast  cancer  detection,  the  boundary  conditions  for  prostate  cancer  detection  are 
non-trivial.  The  effect  of  which  is  obvious  in  the  reconstructed  displacement  and  the  strain  images. 
These  non-uniform  boundary  conditions  result  in  non-uniform  and  non-symmetric  internal 
displacements,  which  can  be  seen  in  the  axial  displacement  image.  This  results  in  strain 
concentration  artifacts  along  the  top  most  boundary  of  the  strain  image.  The  non-uniform  boundary 
conditions  also  cause  a  non-linear  decay  in  strain  with  increasing  depths.  The  isolated  saturation  of 
strain  in  the  lower  quarter  of  the  strain  image  is  due  to  the  out-of-plane  motion,  which  cannot  be 
tracked  with  linear  array  transducer.  Two-dimensional  arrays  may  probably  solve  these  problems, 
which  in  turn  may  result  in  substantial  improvement  in  the  quality  of  strain  images.  It  is  also 
important  to  note  that  no  post-processing  was  done  on  the  strain  images.  In  spite  of  the  unique 
challenges  associated  with  the  elasticity  reconstruction  process  of  the  prostate 
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Fig.  5.  The  lesion  detection  process.  The  lesion  is  barely  identifiable  in  the  original  B-mode  image 
(A).  A  displacement  image  is  illustrated  in  (B).  The  lesion  is  clearly  identifiable  in  the  strain  image 
(C). 


The  3D  volume  of  the  detected  lesion  was  calculated  after  segmenting  the  2D  contour  in  each 
slice  and  estimating  the  inter-slice  distance.  A  multiple  slice  view  and  a  3D  surface  view  of  the 
identified  lesion  were  rendered  based  on  the  3D  dataset  as  shown  in  Fig.  6.  The  volume 
measured  in  the  elastography  and  3D  reconstruction  is  approximately  339±11  mm3  less  than  15  % 
volumetric  error  from  the  volume  of  300  ±30  mm3  measured  using  Archimedes  principle. 


Fig.  6.  3D  reconstruction  of  the  lesion.  A:  The  multiple  2D  image  slices  acquired  were  rendered 
in  three  dimensions.  The  lesion  was  shown  darker  than  surrounding  tissue.  B:  The  multiple  slices 
of  elasticity  images  with  the  detected  lesion  shown  in  black.  C:  The  prostate  lesion  detected 
from  the  ultrasound  images.  The  lesion  was  segmented  from  the  surrounding  tissue  and  its  3D 
surface  was  rendered. 


The  results  shown  above  indicate  that  the  finer  resolution  of  the  new  higher  frequency  transducer 
yield  superior  imaging  and  3D  reconstruction  results.  The  new  transducer  operates  at  up  to  14 
MHz.  The  transducer  previously  used  (as  reported  last  year)  used  8  MHz  imaging  and  a  similar 
aperture  dimension.  Taken  together,  our  in  vitro  results  make  us  well  place  to  commence  initial 
patient  studies  once  our  Army  and  Institutional  Review  Board  permission  is  obtained.  The 
paperwork  associated  with  these  studies  is  being  completed  at  the  time  of  writing. 


Image  Processing  and  Quantification 

We  have  progressed  to  a  accurate  3D  surface  rendering  from  2D  slices  by  implementing  a  3D 
gradient  vector  flow  (GVF)  snake  algorithm  [11],  The  2D  and  3D  GVF  snake  algorithm  relies  on  the 
edges  to  be  well  delineated  and  contrast  between  the  various  regions  to  be  well  defined  in  each  2D 
ultrasound  slice  to  produce  a  surface  that  resembles  the  actual  scanned  object.  In  ultrasound, 
images  are  affected  by  a  granular  pattern  commonly  known  as  speckle.  Before  an  accurate  3D 
surface  rendering  can  be  attained  a  preprocessing  despeckling  step  is  needed  to  reduce  the 
variances  in  pixel  values  within  homogeneous  regions  while  contrast  between  distinct  regions  are 
concurrently  enhanced.  We  have  evaluated  a  wide  variety  of  well  known  methods  such  as  the  Nagao 
and  Matsuyama  filter  [19],  the  Lee  filter[20],  the  Frost  et  at.  filter  [21],  the  Kuan  et  al.  filter  [22],  the 
adaptive  weighted  median  filter  proposed  by  Loupas  et  a/.[23],  the  Wiener  filter  [24],  the  SRAD 
proposed  by  Yu  and  Acton  [25],  and  a  novel  stochastically  driven  method  design  specifically  for  3D 
surface  rendering  from  2D  slices  of  the  prostate  and  other  organs.  The  method  we  develop 
specifically  for  the  task  of  this  grant  is  a  stochastically  driven  compression  filter  called  the  squeeze 
box  filter  (SBF).  Our  quantitative  evaluation  using  a  Field  II  [12]  simulated  B  mode  ultrasound  image 
with  contrast  enhancement  performance  determined  by  a  modified  Fisher  discriminant  has  determined 
that  the  newly  developed  SBF  method  outperformed  the  other  methods  and  is  exceptional  in  providing 
the  needed  intra-region  reduction  in  variance  and  inter-region  contrast  enhancement  with 
computational  efficiency.  In  Figs.  7  and  8,  we  show  the  results  of  SBF  and  SRAD  applied  to  a  Field  II 
simulated  image.  It  is  visually  evident  that  the  bright  and  dark  disks  in  the  SBF  result  is  more 
pronounced  than  in  the  SRAD  image.  In  Figs.  9  and  10,  we  show  the  SRAD  and  SBF  results, 
respectively,  using  an  ultrasound  image  of  a  phantom.  Again,  it  is  evident  that  the  results  of  the 
simulation  are  up  held  with  the  results  of  the  actual  ultrasound  phantom  image  in  that  edges  are  well 
preserved  and  contrast  is  better  enhanced  with  the  SBF  than  with  SRAD.  In  Figs.  13,  14,  and  15,  we 
show  the  3D  surface,  the  side  view,  and  the  bottom  view,  respectively,  of  the  rendering  we  attained 
from  a  sequence  of  scans  we  acquired  of  an  egg  phantom.  The  sequence  consists  of  acquiring  a  2D 
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slice  every  millimeter  along  the  long  axis  of  an  egg  phantom.  We  processed  each  slice  with  the  SBF 
despeckling  method.  The  original  unprocessed  middle  slice  is  show  in  Fig.  11.  The  SBF  processed 
middle  slice  is  shown  in  Fig.  12.  The  3D  rendering  was  attain  by  SBF  processing  each  slice  then 
applying  a  3D  GVF  snake  to  attain  the  final  results  shown  in  Fig.  13,  14,  and  15.  It  is  very  evident  that 
our  method  captured  the  essential  size  and  shape  of  the  egg  phantom.  The  volume  estimate  we 
attained  for  the  object  enclosed  by  the  surface  was  only  1 0%  off  of  the  actual  volume  of  the  phantom. 
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Figure  7.  The  result  of  SRAD  on  a  Field  II  simulated  image. 
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Figure  8.  The  result  of  SBF  on  a  Field  II  simulated  image. 
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Figure  10.  The  result  of  SBF  on  an  actual  ultrasound  image  of  a  phantom. 


Figure  1 1 .  Unprocessed  middle  slice  of  the  egg  phantom. 
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Figure  12.  SBF  processed  middle  slice  of  the  egg  phantom 
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Figure  13.  3D  surface  found  by  the  3D  GVF  with  slices  preprocessed  by  SBF. 
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Figure  14.  Side  view  of  the  surface  rendered  by  the  3D  GVF  snake 
with  slices  preprocessed  by  SBF. 


Figure  15.  Bottom  view  of  the  surface  of  the  egg  phantom 
rendered  by  the  3D  GVF  snake  with  slices  preprocessed  by  SBF. 


A  conference  paper  of  our  novel  despeckling  method  has  been  peer  reviewed  and  accepted  for 
presentation  and  publication  at  the  2006  IEEE  International  Symposium  for  Biomedical  Imaging: 
From  Nana  to  Macro  to  be  held  at  Crystal  Gateway  Marriott  (Arlington,  VA  USA)  on  April  6-9,  2006. 
Another  conference  paper  that  focuses  on  improvements  of  the  SBF  by  using  a  dynamically  varying 
window  is  under  peer  review  for  presentation  and  publication  at  the  IEEE  International  Conference 
on  Image  Processing  2006  to  be  held  at  Atlanta  Marriott  Marquis  (Atlanta,  GA  USA)  with  notification 
of  acceptance  due  April  17,  2006.  A  journal  paper  on  the  motivation  and  mathematical  foundations  of 
the  ID  SBF  has  been  submitted  (September  2005)  to  IEEE  Transactions  on  Signal  Processing  and 
is  current  under  going  the  stringent  peer  review  process.  A  journal  paper  on  the  previously  mentioned 
evaluation,  2D  SBF,  and  3D  surface  rendering  is  currently  being  internally  reviewed  and  we  expect  a 
submission  to  an  academically  reputable  journal  in  the  near  future  (February  2006). 

The  future  direction  of  our  research  in  this  regards  will  be  to  increase  the  efficiency  of  our  algorithm, 
improve  the  3D  GVF  snake  to  incorporate  a  priori  know  information  such  as  the  approximate  size  and 
shape  of  the  prostate,  and  compare  our  methods  with  well  established  ground  truth  (we  may  have  to 
establish  this  ground  truth  ourselves). 

KEY  RESEARCH  ACCOMPLISHMENTS 

•  We  have  largely  completed  Aims  1  and  2.  Specifically,  we  have  designed  and  had  fabricated  a  very 
high  resolution  transrectal  ultrasound  transducer  array  for  high  resolution  prostate  imaging. 

•  We  have  integrated  the  new  transducer  with  an  ultrasound  scanner  and  an  automated  injection 
stage  to  realize  an  accurate  elastographic  imaging  device. 

•  We  have  tested  the  3D  and  elastographic  imaging  capability  of  the  transducer  /  scanner. 

•  We  have  made  significant  progress  in  both  speckle  reduction  and  in  prostate  ultrasound  feature 
segmentation  for  improving  and  automating  the  prostate  cancer  diagnostic  process. 

REPORTABLE  OUTCOMES 

Y.  Li,  A.  Patil  and  J.  A.  Hossack,  “High  resolution  three-dimensional  prostate  ultrasound  imaging”. 

Presented  at  SPIE  Medical  Imaging,  San  Diego,  CA,  2006 

Y.  Li,  A.  Patil  and  J.  A.  Hossack,  “Combined  elasticity  and  3D  imaging  of  the  prostate”  Proceedings 
of  2005  IEEE  Ultrasonics  Symposium,  pp. 1435-1438,  2005 

P.  C.  Tay,  S.  T.  Acton  and  J.  A.  Hossack,  “A  Stochastic  Approach  to  Ultrasound 
Despeckling”,  Accepted  for  presentation  at  2006  IEEE  International  Symposium  for 
Biomedical  Imaging:  From  Nana  to  Macro  Arlington,  VA,  2006 

P.  C.  Tay,  S.  T.  Acton  and  J.  A.  Hossack,  “Ultrasound  Despeckling  Using  An  Adaptive 
Window  Stochastic  Approach”,  Submitted  for  presentation  at  IEEE  International  Conference 
on  Image  Processing  2006 

(We  have  journal  papers  planned  and  early  preparation  in  the  areas  of  elastographic  imaging 
and  in  image  speckle  reduction.) 

CONCLUSIONS 

Our  prostate  imaging  approach  combines  using  an  I-Beam  transducer  with  3D  capability, 
elasticity  imaging  and  test  on  a  prototype  using  a  prostate  tissue-mimicking  phantom.  The 
prostate  strain  imaging  performed  here  using  a  slightly  inflated  sheath  over  the  transrectal 
transducer  significantly  enhanced  tumor  visibility  (a  hard  inclusion  in  the  phantom).  (The  lesion 
was  nearly  invisible  in  the  regular  B-mode  image.)  The  I-Beam  transducer  enabled 
reconstruction  of  discrete  2D  image  acquisitions  into  regular  3D  grid  space,  and  thus  the 
tumor  was  rendered  in  3D.  The  volume  calculated  for  this  tumor  had  an  error  of  approximately 
11%  compared  to  the  actual  (independently  determined)  volume.  Additionally,  we  have  made 
significant  progress  in  the  area  of  image  pre-processing  (i.e.  speckle  reduction)  and  in  image 
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quantification  (i.e.  feature  segmentation).  These  image  processing  contributions  significantly 
enhance  the  practical  utility  of  our  technique  since  they  hold  the  promise  of  accelerating  the 
prostate  cancer  diagnostic  task  and  reducing  intra  and  inter  operator  variability.  Reducing 
variability  has  significance  since  serial  analysis  of  cancer  growth  or  remission  is  dependent  on 
accurate  and  repeatable  measures  of  prostate  volume.  Since  we  are  able  to  measure  volumes 
directly,  rather  than  extrapolating  volume  from  a  length  dimension  or  cross-sectional  area,  our 
image  contributions  are  well-matched  and  complement  our  contributions  in  3D  and 
elastographic  imaging. 
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ABSTRACT 

This  work  reports  on  the  application  of  ultrasound  elastography  to  prostate  cancer  detection  using  a  high  resolution 
three-dimensional  (3D)  ultrasound  imaging  system.  The  imaging  was  performed  at  a  relatively  high  frequency  (14 
MHz),  yielding  very  fine  resolution  that  is  optimal  for  prostate  ultrasound  imaging.  The  fine  resolution  achieved  aids 
in  locating  smaller  lesions  than  are  normally  detectable.  Elasticity  was  measured  with  a  quantitative  and  automatically 
controlled  “Synthetic  Digital  Rectal  Examination  (SDRE)”  wherein  a  smoothly  increasing  force  was  applied  by 
injecting  water,  controlled  by  an  electronic  syringe  pump,  into  a  latex  cover  over  the  transrectal  transducer.  The  lesion 
identified  as  stiffened  tissue  was  visually  enhanced  by  colorizing  and  superimposing  it  over  the  conventional  B-mode 
image.  Experimental  results  using  a  tissue-mimicking  phantom  demonstrated  that  the  reconstruction  accuracy  of  the  I- 
Beam  transducer  resulted  in  less  than  15%  volumetric  error.  Thus,  this  high  resolution  3D  prostate  elastography  is 
possible  and  may  provide  reliable  and  accurate  determination  of  the  size  and  the  location  of  cancers,  which  may  result 
in  improved  specificity  and  sensitivity  of  cancer  detection. 

Keywords:  prostate  cancer,  elastography,  3D  ultrasound  imaging 


1.  INTRODUCTION 

Prostate  cancer  is  the  second  most  prevalent  malignant  cancer  in  North  American  men.  In  2004,  approximately  230,110 
new  cases  of  prostate  cancer  occurred  in  the  US.  According  to  the  American  Cancer  Society  statistics,  29,900  prostate 
cancer  related  deaths  were  anticipated  in  2004  [2].  Systematic  multi-core  biopsy  of  the  prostate,  consisting  of  the 
acquisition  of  ten  or  more  biopsy  cores  distributed  throughout  the  prostate,  is  the  best  current  test  for  detecting  prostate 
tumors,  but  still  exhibits  a  sensitivity  significantly  less  than  ideal  [3].  Most  commonly  used  approaches  for  prostate 
cancer  include  the  digital  rectal  examination  (DRE),  and  the  prostate  specific  antigen  (PSA)  blood  test.  However,  the 
usefulness  of  DRE  is  fundamentally  limited  by  its  subjective  nature,  and  the  PSA  presents  a  significant  “gray”  area  (PSA 
concentration  is  between  4-10  ng/ml)  due  to  its  low  sensitivity  and  specificity  [4].  Recent  examples  of  progress  in  the 
field  of  prostate  cancer  detection  include,  “mechanical  imaging”  of  the  prostate  [5]  and  contrast-enhanced  ultrasound  [6, 
7],  but  each  of  these  have  shortcomings  that  limit  their  application  in  clinical  use.  Mechanical  imaging  has  shown 
noticeable  improvement  in  accuracy  and  sensitivity  for  in  vitro  cases,  while  its  use  of  strain  gauges  and  magnetic 
positioner  device  may  become  cumbersome,  or  even  unfeasible,  for  in  vivo  cases.  The  cost  of  contrast  agents  for  prostate 
cancer  diagnosis  accounts  for  almost  half  of  the  cost  of  a  conventional  ultrasound  examination,  and  therefore  represents 
a  considerable  financial  impediment  to  the  widespread  acceptance  of  contrast  in  routine  prostate  cancer  ultrasound-based 
diagnosis.  In  the  recent  years,  elasticity  imaging  has  garnered  attention  as  a  technique  that  reveals  the  tissue  hardness 
and  thus  provides  a  means  that  complements  anatomic  B-Mode  imaging  with  a  map  of  localized  regions  of  abnormal 
stiffness.  Significantly,  cancers  are  frequently  associated  with  local  changes  in  the  tissue  mechanical  properties,  or  tissue 
hardness  [8].  Techniques  like  sonoelasticity  [9]  and  transient  elastography  involve  application  of  low  frequency  shear 
waves  to  a  tissue.  Difference  in  tissue  elasticity  causes  change  in  the  velocity  of  the  shear  waves  through  the  tissue.  This 
information  can  be  used  to  reconstruct  tissue  elasticity  modulus.  Another  ultrasound  based  elasticity  imaging  is 
“elastography”  [10].  Elastography  has  been  extensively  studied  for  cancer  detection,  and  has  produced  promising  results 
both  in  vitro  and  in  vivo  [11]  [12-15].  Elastography  involves  application  of  quasi  static  compression  to  the  tissue;  the 
resultant  tissue  deformation  is  obtained  by  tracking  the  pre-  and  the  post-compression  echo  RF  data.  The  tissue 
deformation,  or  strain,  is  an  inverse  function  of  tissue  elasticity  and  reveals  the  mechanical  properties  of  the  tissue.  This 
information  can  be  used  to  differentiate  cancers  from  non-cancers  at  a  relatively  early  stage  as  the  spatial  resolution  of 


the  technique  is  comparable  with  that  of  B-mode  ultrasound  [16,  17]. 

One  of  the  challenges  in  elastography  is  to  classify  a  detected  cancer  as  malignant  or  benign,  this  may  potentially  help  in 
precluding  inessential  biopsies.  Garra  et  al.  [13]  proposed  a  cancer  classification  approach  based  on  the  size  of  the  cancer 
as  estimated  from  the  elastographic  images.  It  is  well  known  that  the  benign  cancers  manifest  different  shapes  from 
those  of  the  malignant  cancers.  Hence,  shape  and  size  estimation  of  the  cancers  based  on  3D  volume  (tomography) 
imaging  holds  a  key  to  increasing  the  sensitivity  and  the  specificity  of  elastography  (strain  imaging).  No  3D  prostate 
elastography  studies  have  been  reported  to  date.  We  previously  performed  3D  reconstruction  of  the  in-vitro  prostate 
inclusion  using  I-beam  transducer  operated  at  8  MHz  [18].  This  paper  continues  the  in  vitro  work  on  3D  prostate 
elastography  using  an  I-Beam  transducer  at  a  high  frequency  up  to  14  MHz,  providing  higher  image  resolution  and  thus 
more  accurate  measurement  of  tumor  size  and  location,  especially  for  the  low-volume  tumors. 


2.  METHODS 


2.1  Prostate  Phantom 

A  purpose-built  prostate  phantom  was  designed  using  the  method  described  in  Negron  et  al.  [19].  The  phantom  (Fig.  1) 
comprises  of  four  components  (labeled  1  to  4).  A  pea-sized  (approximately  0.3  ml)  lesion  mimicking  inclusion 
(referred  to  as  “lesion”  hereafter),  an  egg-shaped  tissue  mimicking  prostate  (referred  to  as  “tissue”  hereafter)  a 
cylindrical  cavity  (20  mm  in  diameter)  on  top  of  the  egg-shaped  surrounding  to  mimic  the  anal  opening,  where  the 
transrectal  transducer  is  placed  during  screening;  and  an  exterior  gel  component  that  encases  all  the  other  three 
components.  The  lesion  is  made  of  17%  (by  weight)  acrylamide  gel;  the  tissue  and  exterior  component  are  made  of  5% 
acrylamide  gel.  Thus,  the  lesion  is  perceptibly  stiffer  (10  times)  than  both  the  egg  shaped  tissue  region  and  the  exterior 
gel  component.  The  lesion  and  tissue  are  made  echogenic  by  mixing  Sephadex  (Amersham  Pharmacia  Biotech, 
Piscataway,  NJ),  which  generates  speckle  in  ultrasound  imaging.  A  similar  concentration  (by  volume)  of  Sephadex 
was  added  to  both  the  lesion  and  the  tissue,  resulting  in  similar  ultrasound  image  intensity  in  these  two  structures  (The 
lesion  was  made  slightly  brighter  than  the  tissue  in  order  to  assist  navigation  during  scanning.)  A  B-mode  ultrasound 
image  of  the  prostate  is  shown  in  Fig.  1.  The  lesion  in  the  image  is  almost  isoechoic. 


Legend: 

1  -  lesion  mimicking 
inclusion 

2  -  prostate  tissue 
mimicking 
surrounding 

3  -  rectal  cavity 

4  -  exterior  casing. 


Fig.  1 .  Left:  Schematic  of  the  transrectal  design  of  prostate  phantom.  Right:  An  ultrasound  image  of  the  phantom. 


2.2  Imaging 

A  Siemens  Sequoia  512  scanner  (Siemens  Medical  Solutions,  Mountain  View,  CA)  and  a  linear  phased-array  transrectal 
transducer  operating  at  14  MHz  were  used  in  this  work.  Figures  2  and  3  show  a  photograph  of  the  new,  dedicated,  high 
resolution,  and  3D  capable  “I-Beam”  [20]  transducer  that  was  specifically  designed  for  and  used  in  this  work.  The  I- 
Beam  is  a  modified  ID  linear  array  transducer,  which  allows  for  simultaneous  fast  acquisition  and  tracking  of  the  2D 
planar  images  as  the  transducer  sweeps  in  the  elevational  direction.  This  transducer  comprises  192,  8-14  MHz  elements 


on  a  0.2  mm  pitch  for  imaging  and  two  pairs  of  32  elements  on  a  0.2  mm  pitch,  for  tracking  in  3D  space.  The  I-Beam 
approach  has  been  previously  demonstrated  to  produce  two  standard  deviation  accuracy  in  the  reconstructed  (i.e. 
elevational)  direction  of  4.6%  [20].  The  high  frequency  and  long  aperture  (39.2  mm)  of  this  transducer  enables  possibly 
the  highest  prostate  oriented  ultrasound  imaging  resolution  yet  achieved.  We  anticipate  that  this  high  resolution 
transducer  is  capable  of  more  accurately  measuring  the  size  and  location  of  the  small  lesions  and  we  plan  to  test  it  in  vivo 
in  the  near  future  subject  to  the  requisite  institutional  human  subjects  protection.  In  addition,  the  high  image  resolution 
also  offers  a  potential  for  capturing  lesion  volumes  that  go  undetected  in  low  resolution  imaging  systems.  The  I-Beam 
transducer  is  also  uniquely  matched  to  this  particular  application  for  the  following  reasons:  a)  the  trans-rectal  probe  with 
the  tracking  mechanism  near  the  transducer  prevents  any  numerical  ill-conditioning,  b)  The  I-Beam  transducer  estimates 
the  relative  tissue  motion  rather  than  absolute  tissue  motion,  which  facilitates  efficient  and  accurate  measurements  even 
if  there  is  patient  motion  of  the  type  that  can  defeat  a  3D  system  that  uses  a  fixed  origin  for  3D  positioning  (eg. 
Magnetic-based,  optical-based,  or  articulated  arm-based  positioner.) 

A  quantitative  and  controllable  “Synthetic  Digital  Rectal  Examination”  approach  based  on  automatic  water  inflation  was 
used  to  generate  a  quasistatic  stress  and  to  produce  tissue  deformation.  The  syringe  volume  of  60  ml  was  chosen  to 
provide  sufficient  water  to  compress  and  deform  the  rectal  wall  thus  providing  an  optimal  tissue  strain.  The  injected 
water  induced  uniform  stress  through  the  depth  of  the  tissue,  which  minimized  the  artifacts  due  to  non-uniformities  in 
the  applied  stress.  Other  groups  have  also  adopted  similar  balloon  inflation  method,  however,  those  were  operated 
manually  [  1 4] [  1 5] .  Briefly,  water  was  injected  with  a  syringe  and  ducted  via  a  small  tube  to  the  active  surface  of  the 
transducer  that  was  covered  by  a  latex  condom  (Instead  of  securing  a  tube  from  outside,  an  injection  port  was  configured 
into  the  transducer.)  A  programmable  syringe  pump  (PHD  2000,  Harvard  Apparatus,  Holliston  MA)  was  used  to 
motorize  the  syringe  thus  minimizing  jitter  due  to  undesired  manual  motion. 


Fig.  2.  The  transrectal  transducer  is  covered  with  a  latex  condom.  Water  was  inflated  by  the  syringe  during  imaging 


Fig.  3.  I-Beam  transducer  -  possessing  a  main  imaging  array  in  the  center  and  a  tracking  array  in  each  end 


Six  elevational  slices  each  comprising  of  100  IQ  data  frames  were  acquired  during  the  in-vitro  experiment,  such  that  the 
incremental  strain  between  the  consecutive  frames  was  0.04  %.  Total  applied  strain  between  the  first  and  the  last  frame  was 
approximately  4  %.  An  example  B-mode  image  obtained  by  the  transrectal  I-Beam  transducer  is  shown  in  Fig.  4. 


Fig.  4.  A  B-mode  image  acquired  with  the  I-Beam  transrectal  transducer  from  the  prostate-mimicking  phantom, 
showing  the  layout  of  image  planes  formed  by  ‘Tracking’  arrays  and  ‘Imaging’  array.  ‘T’  -  tracking  array,  ‘I’  - 
main  imaging  array.  Delineated  is  a  lesion-mimicking  inclusion. 


2.3  Elasticity  Calculation  and  3D  Reconstruction 

An  elasticity  map  was  calculated  using  the  acquired  data.  The  data  were  filtered  using  a  low  pass  filter  to  eliminate  jitter 
and  electronic  noise.  Six  pairs  of  frames  with  differential  strain  of  2%  were  tracked  using  a  time-domain  cross¬ 
correlation  technique  on  complex  IQ  data.  Sub  sample  precision  was  acquired  in  the  delay  estimates  by  using  a  quadratic 
fit  to  the  cross-correlation  function.  A  search  window  of  approximately  5  wavelengths  was  used  for  time-delay 
estimation.  Lateral  motion  was  tracked  using  the  technique  described  by  Lubinski  et  al.  [21].  Displacement  estimated 
from  these  six  different  renditions  were  averaged  to  eliminate  uncorrelated  random  noise.  Strain  estimates  were  obtained 
by  taking  the  local  gradient  of  the  displacement  image  (in  axial  direction).  The  above  process  was  repeated  by  sweeping 
the  transducer  in  the  elevational  direction  and  six  elevational  slices  were  obtained.  In  the  elevation  direction,  the  inter¬ 
slice  distance  was  estimated  with  a  block  matching  approach  based  upon  the  minimum  sum  of  absolute  differences 
(MS AD)  algorithm  of  the  I-Beam  ‘Tracking’  data.  The  inter-slice  distance  was  found  to  be  slightly  larger  in  the  deeper 
portion  than  in  the  shallower  portion,  indicating  that  the  transducer  was  rotated  by  small  angles  between  measurements. 
In  the  measurements,  the  slices  were  sampled  with  an  interval  no  greater  than  2  mm,  and  the  rotation  angle  increments 
were  less  than  2  degrees  in  order  to  ensure  that  the  elevational  motion  could  be  accurately  tracked.  Once  the  elevational 
motion  of  these  blocks  was  calculated,  the  acquired  image  slices  were  interpolated  on  to  a  regular  3D  grid  in  Matlab 
enabling  3D  volumes  to  be  rendered. 


3.  RESULTS 

Figure  5  illustrates  the  lesion  detection  process.  The  lesion  was  detected  based  on  the  amplitude  of  the  strain  image.  A 
threshold  technique,  after  a  5x5  low  pass  Gaussian  filtering  to  decrease  noise,  was  applied  on  the  strain  image  using  a 
threshold  value  of  0.8%.  Low  strain  regions  (i.e.  resulting  from  high  stiffness)  were  rendered  using  a  translucent  green 
mask,  superimposed  onto  the  B-mode  image  so  as  to  maintain  the  higher  spatial  resolution  of  the  original  B-Mode  data 
while  simultaneously  highlighting  elastically  anomalous  tissue. 
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Fig.  5.  The  lesion  detection  process.  The  lesion  is  barely  identifiable  in  the  original  B-mode  image  (A).  A  displacement  image  is 
illustrated  in  (B).  The  lesion  is  clearly  identifiable  in  the  strain  image  (C).  Figure  (D)  is  a  hybrid  B-mode  strain  image  in  which  the 
B-Mode  image  is  overlaid  with  a  binary  translucent  mask  outlining  a  region  of  enhanced  stiffness  -  as  determined  by  low  strain. 

The  3D  volume  of  the  detected  lesion  was  calculated  after  segmenting  the  2D  contour  in  each  slice  and  estimating  the 
inter-slice  distance.  A  multiple  slice  view  and  a  3D  surface  view  of  the  identified  lesion  were  rendered  based  on  the 
3D  dataset  as  shown  in  Fig.  6.  The  volume  measured  in  the  elastography  and  3D  reconstruction  is  approximately 
339  +  11  mm3  less  than  15%  volumetric  error  from  the  volume  of  300  ±30  mm3  measured  using  Archimedes  principle. 


15 


c 


Fig.  6.  3D  reconstruction  of  the  lesion.  A:  The  multiple  2D  image  slices  acquired  were  rendered  in  three  dimensions.  The  lesion 
was  shown  darker  than  surrounding  tissue.  B:  The  multiple  slices  of  elasticity  images  with  the  detected  lesion  shown  in  black.  C: 
The  prostate  lesion  detected  from  the  ultrasound  images.  The  lesion  was  segmented  from  the  surrounding  tissue  and  its  3D  surface 
was  rendered. 


4.  DISCUSSION 

The  approach  of  water  inflation  driven  by  a  quantitative  and  automatic  syringe  pump  proposed  in  this  paper  offers  the 
benefits  of  uniform  stress  and  reduces  the  non-uniform  stress  artifacts  and  out-of-plane  problem  due  to  manual  motion. 
In  addition,  it  also  provides  an  opportunity  to  improve  elastographic  signal-to-noise  ratio  and  contrast-to-noise  ration 
by  averaging  multiple  strain  estimates  from  an  image  sequence,  without  any  loss  of  axial  resolution  [17].  Additionally, 
this  quantitative  process  possesses  a  high  degree  of  repeatability,  which  helps  to  decrease  inter-observer  variation  and 
is  required  for  longitudinal  studies  in  the  setting  of  cancerous  tissue  progressing  or  recessing. 


5.  CONCLUSIONS 

We  proposed  a  prostate  cancer  detection  approach  combining  elasticity  and  3D  imaging,  using  an  ultrasound  system 
with  the  exceptionally  high  precision  complex  data,  and  an  I-Beam  transducer  with  3D  capability.  The  approach  was 
tested  on  a  prototype  using  a  prostate  tissue-mimicking  phantom.  The  prostate  strain  imaging  performed  here  using  a 
slightly  inflated  sheath  over  the  transrectal  transducer  significantly  enhanced  tumor  visibility  (a  hard  inclusion  in  the 
phantom,  the  lesion  was  nearly  invisible  in  the  regular  B-mode  image.)  The  I-Beam  transducer  enabled  reconstruction 
of  discrete  2D  image  acquisitions  (including  sets  of  2D  elasticity  acquisitions)  in  to  3D  grid  space,  and  thus  allowed 
the  simulated  tumor  to  be  rendered  in  3D  despite  the  lack  of  contrast  in  the  original  B-Mode  images. 
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Abstract —  A  3D  volumetric  elasticity  reconstruction  method  is 
proposed  for  prostate  elastography.  The  method  uses 
crosscorrelation  based  I/Q  data  tracking  technique  in  2D  with  I- 
Beam  transducer  to  reconstruct  3D  volumetric  elastograms.  The 
elevational  motion  is  tracked  using  block  matching  based  on 
minimum  of  the  sum-of-absolute  differences  between  the 
successive  elevational  frames.  3D  elasticity  reconstruction 
estimates  the  3D  shape  and  size  of  cancers,  which  may  be  an 
important  step  towards  their  classification  as  a  malignant  or  a 
benign. 

Keywords-elastography,  elasticity,  prostate,  3D  reconstruction,  I- 
Beam. 


I.  Introduction 

Adenocarcinoma  of  prostate  is  the  most  prevalent 
malignant  cancer  and  the  second  cause  of  cancer  specific 
deaths  in  men.  Each  year  there  are  approximated  220,000  new 
cases  and  28,900  deaths  [1],  Current  best  approaches  for 
detecting  a  prostate  cancer  include  the  digital  rectal 
examination  (DRE),  and  the  prostate  specific  antigen  (PSA) 
blood  test.  PSA  provides  the  clinician  with  an  indication  of  the 
presence  of  a  cancer  and  hints  on  its  overall  development. 
However,  both  PSA  and  DRE  exhibit  low  levels  of  specificity 
and  sensitivity  with  latter  being  a  subjective  approach. 
Additionally,  prostate  cancer  is  known  to  possess  echogenicity 
that  is  similar  to  that  of  the  surrounding  tissue  and  hence  often 
passes  unnoticed  in  a  trans-rectal  ultrasound  scan  [4].  However, 
cancers  are  known  to  be  stiffer  than  the  surrounding  tissue  [11] 
and  hence  can  be  detected  by  imaging  their  elastic  properties  as 
in  ultrasound  elastography  [15].  Ultrasound  elastography 
involves  application  of  an  external  compression  or  stress 
followed  by  tracking  of  the  pre-  and  post-compression  radio 
frequency  (RF)/  In-phase/Quadrature  (I/Q)  echo  data  to 
produce  maps  of  internal  tissue  displacement,  spatial  derivative 
of  which  produces  a  map  of  local  tissue  strains.  The 
advantages  of  this  quantitative  technique  include  a  higher 
sensitivity,  repeatability  for  serial  analysis,  and  a  spatial 
resolution,  which  is  almost  in  par  with  the  conventional 
ultrasound  B-mode  images  [3-6].  Previous  work  on  2D  prostate 
elastography  was  also  reported  by  Souchon  et  al.  [15]  and 
Alam  et  al.  [2].  In  our  work  we  describe  an  in-vitro  3D 
volumetric  prostate  cancer  detection  technique  based  on 
ultrasound  elastography.  The  3D  portion  of  this  work  was 
enabled  using  an  ‘I-Beam’  transducer  [7].  Our  preferred  strain 
display  approach  is  to  colorize  the  underlying  B-Mode  image 
data  according  to  measured  strain  consequently  maintaining  a 
higher  spatial  resolution  of  the  original  B-Mode  data. 


II.  METHODS 

A.  Prostate  Phantom 

The  experiments  were  performed  on  a  laboratory-built 
acrylamide  based  tissue  mimicking  phantom,  prepared  with  a 
protocol  based  on  the  method  by  Negron  [8].  The  phantom 
comprises  four  components  as  shown  in  figure  1.  A  pea-sized 
inclusion  is  made  of  17%  (by  weight)  acryl  amide  gel  to 
simulate  a  hard  lesion.  An  egg-shaped  surrounding  is  made  of 
5%  acrylamide  gel  to  mimic  the  soft  prostate  tissue.  A  round 
cavity  (diameter  20  mm)  on  top  of  the  egg-shaped  surrounding 
holds  the  transrectal  transducer  and  an  outer  body,  also  made  of 
5%  acrylamide  gel,  encloses  all  the  other  three  structures. 
Sephadex  was  added  to  provide  speckle  to  the  surrounding 
tissue  and  the  inclusion,  but  not  to  the  outside  body.  In  the 
figure  1,  on  the  right,  the  cavity  is  obscured  in  the  image 
because  the  soft  acrylamide -based  outer  body,  which  is  not 
echogenic. 


Figure  1:  The  prostate  phantom,  1-cavity  to  put  probe,  2- 
pea  sized  inclusion,  17%  acrylamidel,  3-prostate  tissue,  5% 
acryl  amide,  4-outer  body,  5%  acryl  amide. 

B.  Imaging 

The  A  latex  condom  over  the  transrectal  transducer  was 
sealed  in  a  conventional  manner  on  the  outside  surface  using 
elastic  bands.  The  condom  was  controllably  inflated  with  water 
using  a  syringe.  The  syringe  volume  of  60mL  was  chosen  to 
provide  sufficient  displacement  of  the  rectal  wall  for  a 
measurable  tissue  strain  signal  while  intrinsically  avoiding  the 
risk  of  over  inflation  of  the  condom.  The  pressure  and 
displacement  are  both  well  within  safe  limits  for  in  vivo  use. 
During  the  syringe  inflation,  a  programmable  syringe  pump 
(PHD  2000,  Harvard  Apparatus,  Holliston  MA)  was  attached 


to  the  syringe  to  control  the  amount  and  rate  of  water  injection. 
This  permitted  a  continuous,  and  uniform  inflation  process  that 
enabled  the  image  data  set  to  be  collected  as  one  image 
sequence  capture  process.  A  Siemens  Sequoia  512  scanner 
(Siemens  Medical  Solutions,  Mountain  View,  CA)  was  used  in 
this  study.  Multiple  demodulated  (In-phase/Quadrature,  I/Q) 
radio  frequency  beam-formed  lines  of  acoustic  data  were 
acquired  from  the  ultrasound  scanner  using  a  research  interface 
employing  an  IQ  data  capture  board.  The  data  were  then 
analyzed  offline  on  a  PC. 
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the  higher  spatial  resolution  of  the  original  B-Mode  data  while 
simultaneously  highlighting  elastically  anomalous  tissue.  The 
volumes  calculated  from  the  three  independent  ultrasound 
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Figure  3:  Observed  B-mode  image  and  the  internal  tissue 
displacement. 


Figure  2:  Schematic  setup  for  prostate  elastography. 

The  I-Beam  transducer  allows  for  simultaneous  acquisition 
and  tracking  of  the  2D  planar  images  as  the  transducer  sweeps 
across  the  phantom;  it  is  compatible  with  our  ultrasound 
biplane  trans-rectal  probes.  The  apparatus  allows  high  accuracy 
in  the  dimensional  measurements  with  standard  deviation  of 
around  4.6%  [7].  The  I-Beam  transducer  is  also  uniquely 
matched  to  this  particular  application  for  the  following  reasons: 
a)  the  trans-rectal  probe  with  the  tracking  mechanism  near  the 
transducer  prevents  any  numerical  ill-conditioning,  b)  The  I- 
Beam  transducer  estimates  the  relative  tissue  motion  rather 
than  absolute  tissue  motion,  for  more  efficient  and  accurate 
measurements  [12]. 

Injection  of  water  in  the  latex  condom  results  in  an 
application  of  a  controlled  stress.  This  external  stress  deforms 
the  phantom.  The  resulting  internal  displacements  are  estimated 
by  tracking  the  pre-  and  the  post-compression  I/Q  echo  data 
using  time-domain  crosscorrelation.  The  shift  in  the  peak  of  the 
crosscorrelation  function  corresponds  to  the  local  tissue 
displacements.  The  reconstructed  tissue  displacement  image  is 
differentiated  to  estimate  the  local  strain.  Elevational  motion  is 
tracked  using  an  I-Beam  transducer  and  block  matching  using 
minimum  of  the  sum-of-absolute  difference  technique.  For 
accurate  elevational  tracking,  the  2D  frames  were  sampled 
every  2  mm  and  with  less  than  3  degrees  of  angular  rotation. 
The  obtained  2D  slices  were  interpolated  thus  enabling  3D 
volume  rendering  using  Matlab  (Natick,  MA). 

III.  RESULTS 

Figures  3  and  4  illustrate  the  lesion  detection  process.  The 
displacement  and  strain  maps  were  calculated  from  the  I/Q  data 
as  discussed  above.  The  reconstructed  elastograms  were 
filtered  using  a  3  X  3  Gaussian  filter  and  time  averaged  using 
five  independent  renditions.  Low  strain  regions  (i.e.  resulting 
from  high  stiffness)  were  rendered  using  a  translucent  green 
mask,  superimposed  onto  the  B-mode  image  so  as  to  maintain 


measurements  of  this  inclusion  were  258,  267  and  274  pi 
corresponding  to  a  mean  value  of  266  pi  with  a  standard 
deviation  of  8  pi.  The  volume  of  the  actual  inclusion  was 
measured  to  be  approximately  300  pi  using  Archimedes’s 
principle  in  a  graduated  cylinder.  The  estimated  standard 
deviation  of  the  reconstructed  inclusion  volume  was  ±30  pi 
(10%). 

strain  %  colorized  lesion 


Figure  4:  Reconstructed  strain  image  with  the  strain 
superimposed  on  the  B-mode  image. 

Nevertheless,  this  early  result  is  encouraging  when  taking 
account  of  the  fact  that  errors  in  all  three  orthogonal 
directions  can  compound  to  degrade  the  final  volumetric 
accuracy.  Furthermore,  in  previous  studies  using  a  transducer 
with  similar  acoustic  characteristics  to  the  one  used  here,  we 
obtained  a  standard  deviation  of  4.6%  in  the  reconstructed 
(transducer  elevation)  dimension  [7].  Figures  5  and  6 
illustrate  the  reconstructed  volume  and  the  original  volume  of 
the  inclusion,  respectively. 


IV.  DISCUSSION 

In  this  work  we  have  used  the  base  band  or  I/Q  data,  which 
decorrelates  much  slower  than  the  RF  data  thereby  enabling  the 
imaging  of  higher  induced  strains  though  at  a  relatively  lower 
resolution  [16].  Also,  3D  volume  rendering  may  help  in  better 
estimating  the  size  and  the  shape  of  the  detected  tumor,  which 
may  be  of  great  value  in  determining  whether  the  tumor  is 
malignant  or  benign  [10].  In  addition,  we  have  acquired  a 
transducer  with  higher  center  frequency,  which  will  improve 


Figure  5:  The  re-construeted  3D  volume  of  the  pea-shaped 
inclusion. 


the  spatial  resolution  of  the  elastographic  strain  images  in  the 
future  studies.  Figure  7  illustrates  the  comparison  of  two  B- 
mode  images  obtained  from  the  current  (8  MHz  ,  spatial 
resolution  of  0.18  mm)  and  the  new  (14  MHz,  spatial 
resolution  of  0.1  mm)  phased-array  transducer.. 


Figure  6:  The  original  dimensions  of  the  pea-shaped 
inclusion. 

It  is  important  to  note  that  the  speckle  pattern  in  the  14  MHz 
transducer  is  finer  than  that  in  the  8  MHz.  Also,  prostate 
elastography  requires  imaging  of  the  shallower  part  of  the 
prostates  and  hence  depth  dependent  attenuation  is  minimized 
with  a  gain  in  spatial  resolution. 


Figure  7:  The  B-mode  images  obtained  from  the  8  MHz 
(left)  and  14  MHz  (right)  transducer. 


V.  CONCLUSIONS 

The  prostate  strain  imaging  performed  here  using  a  slightly 
inflated  sheath  over  the  transrectal  transducer  significantly 
enhanced  tumor  visibility  (a  hard  inclusion  in  the  phantom). 
(The  inclusion  was  nearly  invisible  in  the  regular  B-mode 
image.)  The  I-Beam  transducer  enabled  reconstruction  of 
discrete  2D  image  acquisitions  onto  regular  3D  grid  space,  and 
thus  the  tumor  was  rendered  in  3D.  The  volume  calculated  for 
this  tumor  had  standard  deviation  of  approximately  11% 
compared  to  the  actual  (independently  determined)  volume. 
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Abstract 

Surface  rendering  of  three  dimensional  (3-D)  objects  from  two  dimensional  (2-D)  images  require  features  such 
as  edges  to  be  accurately  delineated  and  contrast  between  differing  regions  to  be  well  pronounce.  In  this  regard, 
images  produced  by  an  ultrasound  system  are  adversely  hampered  by  a  stochastic  process  known  as  speckle.  When 
a  set  of  2-D  ultrasound  slices  are  used  to  autonomously  render  3-D  surfaces,  the  ambiguities  caused  by  the  speckle 
phenomena  prohibit  an  accurate  depiction  of  an  object’s  surface.  In  this  paper,  we  objectively  evaluate  a  variety 
of  known  despeckling  filters  to  be  used  as  a  preprocessing  step  to  remove  or  reduce  speckle  from  each  slice  prior 
to  applying  a  3-D  gradient  vector  flow  active  contour  to  determine  the  surface  of  an  object.  We  provide  a  novel 
efficient  despeckling  method  that  is  well  suited  for  this  pre-rendering  process.  This  novel  despeckling  technique 
visually  displays  excellect  contrast  enhancement  in  both  actual  and  simulated  ultrasound  images.  A  quantitative 
evaluation  on  simulated  ultrasound  images  and  using  a  relative  contrast  performance  enhancement  metric  verifies 
our  qualitative  evaluation.  More  importantly,  an  autonomous  3-D  rendering  using  our  novel  despeckling  method 
yields  an  excellent  approximation  to  the  object’s  actual  surface.  Excellence  in  the  sense  the  size,  shape,  and  volume 
estimate  of  the  object  enclosed  by  the  surface  rendering  accurately  reflect  the  size,  shape,  and  volume  of  the  actual 
object,  resp. 
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I.  Introduction 

An  accurate  3-D  rendering  of  organs  such  as  the  heart,  kidney,  prostate,  liver,  etc.,  would  in  many 
ways  be  a  beneficial  aid  to  health  practioners  in  diagnosing  ailments  or  assessing  the  health  of  that  organ. 
These  3-D  renderings  are  possible  by  layering  2-D  images  or  slices  attained  from  an  incremental  scan  to 
construct  a  3-D  matrix.  A  3-D  surface  of  the  object  can  be  formed  from  the  3-D  matrix,  which  hold  edge 
information  contained  in  the  2-D  slices.  A  straightforward  method  to  determine  the  3-D  surface  from  the 
3-D  matrix  of  2-D  slices  is  to  implement  the  3-D  version  of  the  gradient  vector  flow  (GVF)  active  contour 
(i i.e .  snake).  The  GVF  active  contour  of  Xu  and  Prince  found  in  [1]  is  robust  in  determining  the  surfaces 
of  objects  that  are  in  great  generality  star  convex.  The  GVF  snake  derives  pixelwise  gradient  information 
from  either  the  binary  or  grayscale  edge  map.  A  diffusion  of  the  edge  map  gradients  form  a  vector  flow 
field.  From  some  initial  set  of  snaxels  (short  for  snake  elements  and  are  defined  as  pixels  on  the  current 
contour),  the  GVF  iteratively  forces  the  snaxels  and  the  interpolated  surface  to  converge  towards  strong 
(well  defined)  edge  or  border  locations.  A  surface  interpolation  using  the  final  set  of  snaxels  determines 
the  3-D  surface  of  the  object. 

Before  we  can  apply  the  3-D  GVF  snake  to  the  3-D  matrix  composed  of  2-D  ultrasound  images,  the 
speckle  in  each  ultrasound  image  must  be  reduced  so  that  a  faithful  edge  map  can  be  determined  because 
the  GVF  is  dependent  on  this  edge  map.  Although  the  ultrasound  images  are  affected  by  a  stochastic 
process  seen  as  a  granular  pattern  commonly  known  as  speckle,  the  use  of  ultrasound  as  an  imaging 
modality  is  preferred  for  the  following  reasons: 

.  only  safe  non-ionized  sound  waves  are  used  in  the  scanning  process; 

•  the  portability  of  the  hardware  is  advantageous  for  telemedical  purposes  or  imaging  patients  where 
transporting  them  is  problematic; 

.  cost  is  inexpensive  when  compared  to  other  medical  imaging  modalities;  and 
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.  real-time  or  near  real-time  functional  information  such  as  blood  or  tissue  velocities  can  easily  be 
attained  by  taking  advantage  of  the  Doppler  effect  of  the  sound  waves. 

Goodman  in  [2]  offers  a  physical  explanation  of  laser  speckle  as  due  to  the  roughness  with  respect  to 
the  wavelength  of  the  laser  on  the  surface  being  imaged.  The  cause  of  laser  speckle  can  be  mathematically 
models  as  a  sum  of  a  large  number  of  complex  phasor,  a.k.a.  a  complex  random  walk.  These  complex 
phasors  can  have  a  constructive  or  destructive  relationship  with  each  other.  Thus  bright  and  dark  points  can 
be  observed  in  close  proximity  of  each  other.  The  relationship  and  relevance  of  acoustic  speckle  to  laser 
speckle  is  given  by  Abbott  and  Thurstone  in  [3].  More  relevant  and  detailed  descriptions  of  ultrasound 
image  formation  and  the  statistics  of  ultrasound  speckle  can  be  found  in  [4],  [5].  From  the  experimentation 
performed  in  [4],  the  intensity  or  pre-logarithm  compressed  enveloped  detected  amplitudes  can  be  modeled 
as  the  multiplicative  model  given  in  equation  (1) 

J(n,  m)  =  P(n,  m )  *  (I(n,  m)r}x  ( n ,  m))  (1) 

where  rjx(n,m)  is  statistically  independent  of  I (n.  rn)  or  as  the  additive  model 

J(n,  m)  =  P(n,  m)  *  I (n,  m)  +  r}+(n,  m)  (2) 

where  r]+(n,m )  is  statistically  dependent  to  I(n,m).  In  both  equations  (1)  and  (2),  the  point  spread 
function  (PSF)  of  the  ultrasound  imaging  system  is  denoted  by  the  P(n,  m)  and  maybe  spatially  varying. 

The  smoothing  of  speckle  and  preservation  of  edges  are  in  a  general  sense  opposite  processes,  which 
are  difficult  to  accomplish  concurrently.  This  problem  is  not  only  directly  relevant  to  enhancing  ultrasound 
images,  but  also  relevant  to  enhancing  synthetic  aperture  radar  (SAR),  optical  laser,  and  other  imaging 
methods.  A  wide  variety  of  methods  have  been  proposed  to  address  speckle  removal  or  reduction, 
commonly  referred  to  as  despeckling,  with  success  dependent  on  improvements  to  a  post-despeckling 
application  or  visual  interpretation.  When  addressing  speckle  as  unwanted  noise  non-linear  adaptive 
filtering  [6]— [  1 2]  and  anisotropic  diffusion  [13],  [14]  methods  have  been  proposed.  Although  other 
more  recent  types  of  methods  like  multiscale  thresholding  [15],  Bayesian  multiscale  method  [16]  or 
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improvements  on  existing  methods  such  as  those  in  [17]  have  been  proposed,  we  restrict  our  survey  to 
the  Wiener  filter  [18];  the  non-linear  adaptive  methods  of  Nagao  and  Matsuyama  [6],  the  Lee  filter  [7], 
the  Frost  et  al.  filter  [8],  the  Kuan  el  al.  filter  [10];  and  the  anisotropic  diffusion  method  SRAD  [14]. 
Later  in  this  paper,  a  quantitative  evaluation  on  the  performance  of  these  filters  in  improving  the  contrast 
enhancement  will  be  assessed. 

The  problem  of  attaining  the  noise  free  unblurred  image  I  from  a  detected  image  J,  modeled  as  in 
equations  (1),  (2),  or  others,  is  generally  considered  ill  posed,  that  is,  the  solution  may  not  be  unique  or 
impossible  to  verify.  To  state  a  well  posed  problem,  either  the  restoration  is  to  be  accomplished  under  some 
constraints  on  I,  i.e.  regularization  [19],  [20],  or  the  ideal  unblurred  noise  free  image  I  is  known  a  priori. 
The  experiments  and  results  using  simulated  ultrasound  images  given  later  in  this  paper  will  be  in  the  spirit 
of  the  latter.  We  will  evaluate  the  performance  of  various  sharpening  and/or  denoising  algorithms  based 
on  the  a  priori  knowledge  of  the  desired  image.  The  general  goal  of  an  image  restoration  or  segmentation 
algorithm  should  be  to  decrease  the  variance  in  a  homogeneous  region  while  distinct  regions  should  be 
well  defined.  A  novel  metric,  which  is  a  generalization  of  the  Fisher  discriminant,  will  be  given  to  quantify 
how  well  various  algorithms  achieve  this  general  goal.  An  evaluation  of  various  despeckling  algorithms 
using  simulated  ultrasound  images  will  be  presented.  This  evaluation  will  provide  support  for  the  use  of  a 
novel  stochastically  driven  compression  filter  as  a  preprocessing  step  prior  to  3-D  surface  rendering  using 
a  3-D  GVF  snake.  We  substantiate  our  results  by  comparing  the  3-D  surface  rendering  produced  by  the 
3-D  GVF  snake  using  a  human  optimized  SRAD  followed  by  a  morphological  hole  filling  operation  with 
the  surface  rendering  produced  by  the  3-D  GVF  snake  using  our  proposed  efficient  compression  method. 

II.  Various  Despeckling  Filters 

A.  Wiener  Filter 

Suppose  the  detected  image  is  given  as  the  additive  model  in  equation  (2)  where  the  PSF  P  is  known 
or  estimated.  Applying  the  Wiener  filter  as  described  in  [18]  to  the  detected  image  J  results  in  an  ideal 
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(or  approximately  ideal)  low  pass  filtered  version  of  the  ideal  image  I  polluted  with  the  noise  component 
rj+(n,m),  provided  that  the  PSF  P  is  an  low  pass  filter.  Wiener  filtered  image  J  is  determined  as 

J  j  =  PJ  =  P  |  PI  +  rj+ 1  (3) 

D  FT  ^  ^  ^  ^  ^ 

where  P  < — >  P  (i.e.  P  and  P  are  discrete  Fourier  transform  pairs),  likewise  for  J,  I,  p+  and  J,  I,  r/+ , 


resp. 


P*  (k,l) 


if  P(k,l)  ±  0 


P(k,l)=< j  lP(M)l  +CT'+ 

1  otherwise, 


(4) 


and  afl+  is  the  variance  of  the  noise.  The  fortuitous  side  of  applying  the  Wiener  filter  is  that  edges  are  in 
a  generally  relative  sense  enhanced.  Unfortunately,  the  Wiener  filter  indiscriminately  enhance  the  speckle 
present  within  a  homogeneous  region. 

It  should  be  noted  that  the  Wiener  filter  is  applied  in  the  DFT  domain.  This  requires  a  priori  knowledge 


of  the  PSF,  which  for  ultrasound  images  is  spatially  varying.  Even  if  the  spatially  varying  PSF  is  known  or 
could  be  accurately  calibrated,  incorporating  a  spatially  varying  PSF  into  the  Wiener  filtering  framework 


is  problematic  and  will  be  topics  for  our  future  research  papers. 


B.  Nagao  Filter 

Nagao  and  Matsuyama  in  [6]  proposed  a  recursive  edge  preserving  smoothing  algorithm.  One  iteration 
of  this  algorithm  replaces  each  pixel  value  with  the  mean  of  some  segment  Wk  originating  from  J(n,m ) 
and  the  variance  of  Wk  is  the  minimum  variance  attained  from  a  set  of  variances  in  all  eight  directions. 
Precisely, 

J(n,m)  =  E(wk),  (5) 

where  E(-)  is  the  expected  value  operator, 

var(wk )  =  min  {var(wi)  \  for  i  =  0, 1,  2, . . . ,  7} 

and  w0,wi,u>2,  ■  ■  •  ,u>7  are  equal  length  segments,  which  originate  from  J(n,m)  and  span  all  eight 
directions.  This  recursion  is  allowed  to  continue  until  convergence  in  most  of  the  pixel  values  is  established. 
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This  method  they  claim  and  empirically  show  in  [6]  enhances  the  edges  of  an  image  in  addition  to 
adaptively  smoothing  the  image. 

Though  this  method  is  simplistic,  it  is  evolutionary  in  which  this  filter  adaptively  uses  local  statistics, 
variances  in  the  eight  different  direction,  to  determine  a  window  to  derive  a  replacement  pixel  value.  In 
addition,  the  iterative  nature  of  this  method  may  offset  some  of  the  problems  with  using  a  fixed  window 
length. 

C.  The  Lee  Filter 

Lee  in  [7]  proposed  methods  to  contrast  enhance  an  image  and  to  restore  an  image  corrupted  by  noise. 
Lee’s  noisy  image  models  in  [7]  are 

J(n,  m)  =  I(n,  m)  +  77(71,  m)  (6) 

for  the  additive  noise  case  and 

J(n,m)  =  I(n,m)rj(n,m )  (7) 

for  multiplicative  noise.  In  addition  to  proposing  a  method  to  contrast  enhance  an  image,  his  paper 
proposes  an  adaptive  filter  to  aggressively  smooth,  via  local  averaging,  in  homogeneous  regions  while 
regions  which  contain  significant  image  features  such  as  edges  are  to  be  left  unmolested. 

To  emphasize  the  significance  of  Lee’s  contribution  in  [7],  the  interrelated  contents  of  his  paper, 
contrast  enhancement,  additive  and  multiplicative  noise  suppression  algorithms  are  briefly  and  thoroughly 
described.  To  enable  contrast  enhancement,  at  each  pixel  Lee’s  algorithm  used  a  linear  rescaling  of  the 
local  mean  summed  with  the  multiplication  of  a  gain  applied  to  the  difference  of  the  pixel  value  with  the 
local  mean  to  determine  the  new  pixel  value.  Formally,  let  J(n,m)  be  the  original  pixel  value  of  some 
image  J,  then  the  new  pixel  value  J(n,m )  is  set  at 

J(n,m )  =  g{jx)  +  k(J(n,m )  —  n),  (8) 

where  /i  is  the  local  mean.  The  function  <?(•)  is  a  linear  rescaling  of  the  mean  that  is  g(fj)  =  a/i  +  b  where 
the  parameters  a,  b  e  M  were  determined  to  allow  the  new  pixel  value  to  utilize  the  full  eight  bit  dynamic 
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range.  As  pointed  out  in  [7],  if  0  <  k  <  1,  then  equation  (8)  determines  a  smoothing  filter,  i.e.,  a  low 
pass  filter.  When  the  gain  k  is  greater  than  one,  then  equation  (8)  attempts  to  “sharpen”  image  features 
that  is  enhance  the  edges. 

If  the  image  is  determined  or  assumed  to  be  polluted  with  zero  mean  additive  white  noise1,  then  the 
gain  k  is  adaptively  chosen  as  a  function  of  the  local  statistics.  To  be  precise,  the  new  pixel  value  J(n,  m) 
is  set  at 

J(n,m)  =  ji  +  k(J(n,m)  —  jj)  (9) 


where  /i  is  the  mean  in  some  window.  The  gain  k  is  determined  as 


k  = 


~  < 

a2 


(10) 


where  a2  is  the  local  variance  in  the  same  window  about  J(n.  m)  that  determined  //  and  a2  is  the  global 


noise  variance.  When  a2  S>  a2  ^  0,  the  gain  parameter  is  less  than  but  approximately  one.  In  which  case 


the  filter  in  equation  (9)  performs  like  an  identity  filter  that  is  J(n,  m)  ~  J(n,  m).  If  the  local  variance  a2 
is  greater  than  but  nearly  equal  to  the  global  noise  variance  a2,  then  J(n,m )  ~  //  and  the  filter  specified 
by  equation  (9)  serves  as  a  low  pass  filter.  Since  the  global  noise  variance  can  only  be  greater  than  or 
equal  to  zero,  the  gain  parameter  k  can  never  exceed  one.  Thus,  in  homogeneous  regions  k  should  be  set 
to  zero  and  equation  (9)  provides  local  smoothing.  When  J(n,m )  is  determined  to  be  an  edge  pixel,  then 
k  should  be  set  to  one  and  the  pixel  is  left  undeteriorated. 

When  the  image  is  determined  or  assumed  to  be  degraded  by  multiplicative  noise  as  in  equation  (7), 
then  Lee  in  [7]  approximates  the  image  as  an  additive  noise  model  of  the  form 


J(n,  m)  =  AI(n,  m)  +  Brj(n,  m )  +  C 


(11) 


where  A,  B,C  e  M  are  chosen  so  that  mean  square  error  between  the  .Jin.  m)  of  the  multiplicative  model 
from  equation  (7)  and  J(n,m )  of  the  additive  model  in  equation  (11)  is  minimized.  Since  the  image 
polluted  by  multiplicative  noise  is  recast  as  an  image  with  additive  noise,  the  adaptive  filter  employed  in 

*Lee  in  [7]  does  not  specify  any  distribution  on  the  noise  values. 
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the  additive  case,  given  in  equation  (9),  can  be  used,  though  the  gain  parameter  k  is  redefined.  Given  an 
a  priori  determination  of  the  mean  and  variance  of  the  stationary  noise,  pv  and  a 2  resp.,  the  adaptive 
gain  parameter  is  determined  as 


k 


p,jQ 

+  Ah2, Q 


(12) 


where  p,j  is  the  mean  of  / (n,  m)  determined  as 


pi  = 


p_ 

AG 


(13) 


within  some  fixed  window  and  p  is  the  local  mean  of  J(n,m )  within  the  same  window.  The  variable  0 
is  determined  as 

rr2  A-  ii. 2 

(14) 


^_cr2  +  ^2  _2 

Q  _9  .  .9 


cr2  +  /r2 

where  cr2  is  the  local  variance  of  J(n,m )  within  some  window. 

Lee’s  modification  of  the  adaptive  approach  of  R.  Wallis2  is  a  significant  contribution  because  his  method 
incorporated  local  statistics  of  an  image  to  determine  the  gain  parameters  of  an  adaptive  filter  to  smooth 
additive  or  multiplicative  noise.  Discouragingly,  the  Lee  filter  relies  on  an  unrealistic  approximation  of  the 
local  mean  of  the  ideal  image  pi  given  in  equation  (13).  This  is  not  favorable,  since  this  approximation 
is  difficult  (if  not  impossible)  to  verify  its  accuracy  and  can  be  considered  a  blind  variable. 


D.  The  Frost  Filter 

Frost,  Stiles,  Shanmugan,  and  Holtzman  in  [8]  derived  that  an  image  recorded  from  a  synthetic  aperture 
radar  system  can  be  modeled  as  a  PSF  applied  to  an  ideal  image  polluted  by  multiplicative  noise.  Their 
image  model  is  given  in  equation  (1).  Their  goal  in  [8]  was  to  apply  an  optimum  minimum  mean  squared 
error  (MMSE)  filter  to  estimate  the  local  regions  of  a  ideal  image  while  avoid  degrading  edge  features.  The 
motivation  of  Frost  et  a/.’s  method  consist  of  determining  a  filter  f(n',  m'),  so  that  in  a  local  homogeneous 
region  of  J(n,m )  and  in  the  presence  of  white  noise,  the  expectation  criteria 

E[(I(n,m)  —  I(n,m))2]  (15) 

“The  reference  to  R.  Wallis’  contributions  is  cited  in  [7], 


is  minimized.  The  term  I(n,m)  is  a  windowed  weighted  sum  about  I  in.  rri).  Formally, 

[A  m 

I  (n,  m)  =  E  E  f(n',  m')I (n  +  n',  m  +  m!)  (16) 

where  [x\  denotes  the  greatest  integer  less  than  or  equal  to  x  and  N'  x  M'  is  the  window  size.  Normally, 
N'  and  M'  are  equal  and  odd.  Under  the  assumption  that  the  noise  is  white  and  the  transfer  function  of 
the  PSF  is  constant  over  some  finite  bandwidth,  they  derived  the  filter  that  minimizes  (15)  as 

f{ri,  rri)  =  /Tae-“|T(n'’m,)l  (17) 


where 


a  —  \  2a  \  —  1  1+  - 


K  is  some  normalizing  constant,  a  is  an  region  dependant  constant,  jir]  is  the  mean  of  the  noise,  <rv  is  the 
standard  deviation  of  the  noise,  /j  is  the  local  mean  in  some  window  about  J(n,m),  and  a  is  the  local 
standard  deviation  of  J(n,m )  within  the  same  window  that  is  used  to  determine  the  mean.  The  window 
size  used  to  determine  the  local  mean  and  local  variance  of  J(n,m)  does  not  necessarily  have  to  equal 


the  weighting  window  size,  N'  x  M' .  The  function  r  :  Z  x  Z 


was  not  explicitly  specified  in  [8]. 


Rather,  a  description  of  r(n,m)  was  simply  given  as  a  monotonically  decreasing  function.  This  can  be 
achieved,  for  example,  when  r(n,  m)  =  \fn 2  +  m?  for  n  =  —  J—L,  and  m  =  — — .  The 
window  dimensions  N'  and  M'  are  usually  odd  and  equal. 


E.  The  Kuan  Filter 

Kuan,  Sawchuk,  Strand,  and  Chavel  in  [9]  develop  a  nonlinear  filter  that  they  claim  in  a  local  neighbor¬ 
hood  provides  a  linear  minimum  mean  square  error  (LLMMSE)  estimate  I  (n,  m )  of  the  ideal  noise  free 
image  I(n,m )  from  the  sensed  or  detected  image  J(n,m).  In  1987,  Kuan  et  al.  published  the  application 
of  the  LLMMSE  estimate  to  image  restoration  in  [10].  The  LLMSE  estimate  of  the  ideal  noise  free  image 
is  given  as 

c o2  —  a2 

- 2 — (J(n,m)  -p.) 


I  (■ n ,  m)  —  n  + 


(19) 
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where  p  is  the  local  variance  in  some  window  about  J(n.  m)  and  u2  is  the  local  variance  in  some  weighted 
window  about  J(n,  m).  The  local  weighted  variance  u2  is  defined  as 

t  iti  m 

^  =  Waf  S  w(n'’  m/)(J (n  +  n'i m  +  m> )  ~  z1)2  (20) 

-'=-m  -'=-m 

where  m! )  =  1,  Ar'  and  M'  are  typically  odd  and  equal. 

n'  ,mf 

F.  The  Adaptive  Weighted  Median  Filter 

Loupas,  McDicken,  and  Allan  in  [11]  develop  a  filter  they  named  the  adaptive  weighted  median  filter 
(AWMF).  The  AWMF  was  specifically  designed  to  suppress  speckle  noise  inherent  to  ultrasound  images. 
It  is  should  be  noted  that  the  motivation  for  the  AWMF  is  motivated  by  the  following  ultrasound  image 
model 

J(n,m )  =  I(n,m )  +  y/l(n,m)r)(n,m).  (21) 


It  is  easily  derived  from  the  image  model  given  in  equation  (21)  that  the  local  variance  of  the  observed 
image  J(n,  m)  in  a  homogeneous  region  is  proportional  to  the  noise  variance,  provided  that  the  //(n,  m)  is 
independent  of  J(n,  m).  Precisely,  let  the  ideal  image  be  equal  to  a  constant  c  in  some  local  neighborhood, 

/ (n,  m)  =  c,  then 

a2  =  cal  (22) 


where  a2  and  a2t  are  the  local  variance  of  the  observed  image  and  the  variance  of  the  wide  sense  stationary 
noise,  resp.  With  this  motivation  Loupas  et  al.  defines  the  AWMF  as 


J(n,  m)  =  median/  J(n  +  n\  m  +  m/),  J(n  +  n,  m  +  m'), 


,  J(n  +  n',  m  +  m!) 


w(n'  ,m') 


(23) 


where  w(n',m')  is  a  nonnegative  integer.  Loupas  et  al.  in  [11]  defined  the  N '  x  N',  where  N '  is  odd, 
window  of  weight  terms  as 


w(n  ,  m!)  =  round . 


nonneg 


»•(().  o)  - 


co2\/n'2  +  m 12 

p 


(24) 
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where  rouridnonneg  { • }  means  round  to  the  nearest  non-negative  integer,  c  is  some  constant,  //  is  the  local 
mean,  and  a2  is  the  local  variance.  The  constant  c  and  the  window  term  w( 0,  0)  determines  the  AWMF’s 
ability  to  preserve  edges. 


G.  Speckle  Reducing  Anisotropic  Diffusion 

Inspired  by  the  time  dependant  heat  diffusion  equation,  Perona  and  Malik  in  [13]  published  a  method  to 
perform  anisotropic  diffusion  on  noisy  images.  In  [13],  the  diffusion  equation  of  an  image  with  continuous 
variables  over  time  J(x,y,t )  is  given  as 

J(x,y,t)  =  div{c(x,y,t)V  J(x,y,t)}  =  c(x,y,t)V2J(x,y,t)  +  Vc(x,y,t )  •  VJ{x,y,t )  (25) 


where  (x,y,t)  G  M3,  div  is  the  divergence  operator,  V  is  the  gradient  operator,  V2  is  the  Laplacian 
operator,  v  •  w  represents  the  dot  product  of  two  vectors,  and  J(x,y,  0)  =  J(x,y).  When  c(x,y,t)  is  a 
constant,  then  equation  (25)  defines  isotropic  diffusion.  For  anisotropic  diffusion,  the  function  known  as 
the  conduction  coefficient  ( a.k.a .  coefficient  of  variation,  diffusion  coefficient)  c(-)  is  given  as 


c{x,  y,  t )  =  exp 


(\\ VJ(x,y, 

1  K 


V 


(26) 


or 


c{x,y,t) 


(\\VJ(x,y,t)\\ 
\  K 


-i 


(27) 


where  K  is  some  constant  and  ||  ■  ||  is  the  norm. 

Yu  and  Acton  in  [14]  proposed  a  anisotropic  diffusion  method  where  the  coefficient  of  variation  is 
determined  as  a  function  of  the  ratio  of  the  instantaneous  variance  to  the  instantaneous  mean  squared.  In 
other  words,  their  speckle  reducing  anisotropic  diffusion  method,  which  they  acronymically  named  SRAD, 
utilizes  the  instantaneous  coefficient  of  variation  (ICOV)  as  a  variable  in  their  conduction  coefficient 
function.  The  ICOV  of  SRAD  determine  whether  a  pixel  should  be  smoothed  or  left  unbludgeoned.  The 
SRAD  algorithm  iteratively  processes  a  nonzero  valued  image  I(x,y )  :=  I(x,y;  0)  according  to 

dl(x  y,  t )  =  div[c{q)VI (x,  y\  t)]  (28) 
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and 


dl(x,y;t) 


dn 


=  0 


(29) 


B{U) 


where  n  is  the  outward  normal  vector  to  /j(Q).  the  border  of  Q.  Equations  (28)  and  (29)  are  known  as  the 
SRAD  partial  differential  equations  (PDEs).  The  diffusion  coefficient  c(-)  is  defined  either  as  the  quotient 

«o(*)  +  ?o(*) 


c(q{x,y;t))  = 


Qo(t)  +  q2{x,y;t) 


(30) 


or  as  the  exponential  function 


claix  «■  t) )  =  exD  ( 1oW  ~ 
[q(X’V’)]  P(  4(t)  +qg(t) 


(31) 


In  both  equations  (30)  and  (31),  if  q(x,y;t )  ~  qo(t),  then  c(q(x,y;t ))  ~  1  and  smoothing  with  respect  to 
equation  (28)  is  enacted.  If  q(x,y;t )  3>  qo(t),  then  the  diffusion  coefficient  is  very  small  and  smoothing 
in  a  local  region  around  (x,  y)  is  averted.  When  at  time  t,  if  (x,  y)  resides  in  a  homogeneous  region,  then 
smoothing  can  be  promoted  by  allowing  q(x,  y,  t )  «  q0(t).  When  (x,  y)  lie  on  an  edge  or  in  the  vicinity  of 
an  edge,  then  defining  q(x,y,t )  3>  qo(t)  would  prohibit  deterioration  of  the  edge.  Yu  and  Acton  defined 
qo(t),  the  coefficient  of  variation  in  fully  developed  speckle,  as  the  ratio  of  the  noise  standard  deviation 
to  the  noise  mean, 


qo(t)  = 


(32) 


and  the  instantaneous  coefficient  of  variation  is  defined  as 


q(x,y;t )  = 


1  ( |v/|  V  _  j_  (v2iy 

2  \  I  I  16  W  / 


(33) 


\  (i  +  d)  (Y))2  ' 

The  standard  deviation  and  the  mean  of  the  noise  in  equation  (32)  are  determined  by  calculating  the  mean 
and  standard  deviation  within  a  homogeneous  region  where  noise  is  prevalent. 

The  coefficient  of  variation  used  in  equations  (26)  and  (27)  rely  solely  on  the  gradient  norm.  The 
inclusion  of  the  Laplacian  operator  in  determining  the  instantaneous  coefficient  of  variation  defined  in 
equation  (33)  and  the  characterization  of  fully  developed  speckle  by  equation  (32)  provides  a  more  robust 
method  to  determine  the  diffusion  coefficient  defined  in  equation  (30)  or  (31)  within  different  possibly 
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disconnected  region.  Lastly,  it  should  be  noted  that  though  we  outlined  the  intuitive  ideas  of  SRAD 
using  continuous  variables  (x,  y,  t )  €  M3,  our  implementation  of  SRAD  used  in  the  evaluation  given  in 
section  IV  of  this  paper  is  the  discrete  form  detailed  in  [14]. 

III.  A  Stochastically  Driven  Method 

With  the  exception  of  the  Wiener  filter,  the  despeckling  methods  described  in  section  II  evaluate  every 
sample  in  a  image  and  adaptively  determines  whether  to  smooth  (locally  average)  or  not.  The  Wiener  filter 
relies  on  a  priori  knowlegde  of  the  PSF  and  the  noise  variance  to  determine  a  convolution  kernel,  which 
is  applied  to  every  pixel  in  the  entire  image.  The  iterative  filtering  method  we  presented  in  this  section 
only  considers  samples  which  are  outliers  of  some  probability  density  function  (PDF)  and  applies  local 
smoothing  to  these  outliers.  The  local  extrema  are  considered  outliers  and  are  not  used  in  the  determination 
of  the  local  mean.  The  choice  of  the  neighborhood  J\f  is  extremely  important,  since  the  mean  of  some 
PDF  is  determined  by  samples  in  J\f .  A  large  neighborhood  fully  contained  in  a  homogeneous  region  will 
produce  a  local  mean  that  is  closer  to  the  mean  of  the  homogeneous  region.  Yet,  a  neighborhood  that  is 
too  large  will  extend  pass  the  homogeneous  region  and  yield  a  local  mean,  which  will  erroneously  reflect 
the  mean  of  the  homogeneous  region.  Each  iteration  of  the  method  currently  being  described  produces  a 
sequence  with  locally  reduced  variance.  The  local  extrema  of  the  new  sequence  are  consider  as  outliers 
and  the  process  is  iterated.  The  steps  of  our  proposed  iterative  method  are  as  follow: 

1)  Each  iteration  i  begins  by  determining  the  set  of  locations  of  local  maxima  and  local  minima.  The 
locations  of  these  extrema  are  defined  by  the  set 

Me  =  {(n,m)  |  J*_i (n,m)  meets  condition  1  or  2  } 

Condition  1:  m)  >  J,_i(n  +  k,  m  +  l ) 

Condition  2:  m)  <  +  k,m  +  l ) 

where  k,l  —  —  1  or  1. 


13 


2)  Without  using  the  local  extrema  values,  our  algorithm  replaces  the  extremum  with  the  local  mean 
taken  from  neighboring  samples.  For  all  (n,  m)  €  Me 

Mn,m)  =  -^jr  Y  Ji~ i(M)  (34) 

where  A/”  is  some  local  neighborhood  of  (■ n,m ),  |A/"|  is  the  cardinality  of  set  A/”,  and  ( n,m )  ^  A/”. 

3)  If  convergence  is  not  attained,  that  is 

^  \Ji_i(n,m)  —  Ji(n,m)\  >  e  (35) 

Vn,m 

for  some  predefined  e  >  0,  then  another  iteration  is  performed.  If  convergence  is  attained,  then  only 
further  trivial  insignificant  improvements  can  be  attained  with  this  filtering  method  and  the  process 
is  stopped. 

By  removing  outliers  at  each  iteration,  this  method  reduces  the  local  variance  at  each  pixel.  In  effect, 
this  method  produces  a  convergent  sequence  of  images  by  squeezing  or  compressing  the  stochastically 
distributed  pixel  values  to  a  limiting  value.  Thus,  we  call  this  stochastically  driven  method  the  squeeze 
box  filter  (SBF). 


IV.  Evaluations  and  Results 

A.  Using  Simulated  Images 

To  evaluate  the  despeckling  performance  of  the  Wiener  filter,  the  various  adaptive  non-linear  filters, 
and  the  anisostropic  diffusion  method  SRAD  against  the  stochastically  driven  SBF  method,  four  templates 
where  created  and  are  denoted  as  T 1 .  T11,  T111,  and  T 1 '  .  These  templates  represent  the  ideal  noise  free 
images  and  are  shown  in  the  top  row  of  Fig.  1.  The  pixel  values  of  the  background  class  of  all  templates 
are  set  at  one.  Templates  T7  and  Tn  consist  of  two  classes:  the  background  class,  shaded  in  black  for 
T7  and  white  for  T ;  1 ;  and  a  class  consisting  of  ten  disks  of  various  size,  shaded  in  white  and  with  pixel 
values  five  for  T7,  and  for  T 1 1  this  class  is  shaded  as  black  and  the  pixel  values  are  set  at  zero.  Templates 
TUI  and  T 1 1  consist  of  three  classes:  the  background  class,  which  is  shaded  in  gray;  a  class  consisting 
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of  five  disks  of  various  diameter  shaded  in  white  with  the  pixel  values  set  at  five  for  TUI  and  zero  for 
T 1  v’ ,  shaded  in  black;  the  third  class  is  composed  of  five  disks  of  various  diameters,  shaded  in  black  and 
with  pixel  values  set  at  zero  for  Tin  and  for  TiV  this  class  is  shaded  in  white  and  with  pixel  values  set 
at  five. 

The  Field  II  simulation  [21]  of  the  four  templates  are  shown  in  the  bottom  row  of  Fig.  1  and  are  denoted 
as  J1,  J11 ,  J111,  and  Jn  .  The  simulations  are  constructed  with  the  transducer  at  the  top  of  the  image.  The 
focus  point  for  each  simulation  is  set  at  70mm  axial  distance  from  the  transducer  and  at  lateral  position 
0mm.  In  each  simulation,  the  spatial  varying  PSFs  are  shown  along  the  column  at  lateral  position  -15mm 
at  axial  distances  40mm,  50mm,  60mm,  70mm,  and  80mm. 

To  evaluate  which  method  in  sections  II  and  III  provide  the  best  improvements  to  an  ultrasound  image 
a  meaningful  quantifiable  measure  and  a  method  to  attain  this  measure  is  needed.  This  measure  should 
indicate  when  different  homogeneous  regions  are  properly  defined.  Additionally,  this  metric  should  account 
for  differences  in  pixel  values  from  the  mean  of  a  class.  A  class  is  taken  to  be  a  collection  of  homogeneous 
regions.  This  measure  in  essence  will  determine  how  well  an  algorithm  despeckles  the  simulated  ultrasound 
image  while  keeping  the  distinct  classes  well  separated. 

The  method  we  propose  to  quantify  the  improvements  made  to  an  ultrasound  image  is  as  follow: 

1)  First,  we  create  the  ideal  image  that  is  a  template  T  with  various  classes  C\ ,  C'2,  C'3, . . . ,  CN. 

2)  A  simulation  of  an  ultrasound  image  J  using  the  template  T  is  accomplished  via  the  Fields  II 
software. 

3)  A  despeckling  algorithm  is  applied  to  the  simulated  image  J.  The  output  image  is  denoted  as  I aig- 

4)  The  means  and  variances  in  each  predefined  classes  of  I aig  are  computed. 

5)  The  ratio  of  the  average  squared  differences  of  the  inter-class  means  to  the  sum  of  the  intra-class 
variances  are  calculated  and  this  quantity  is  denoted  as  Qialg  .  Explicitly,  the  quantity  to  evaluate  the 
performance  of  an  despeckling  algorithm  to  preserve  distinct  classes  and  promote  smoothing  within 
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homogeneous  regions  of  each  class  is  given  as 


E 


Vc, 


q  def  fc/d 

*-alq 


N 


(N-  l)^a 


2 


fc=l 


where 


de/  1 


Vck  =  7777  E  Wn>m)> 

(n,m)GCfc 


|Cfc| 


(36) 


(37) 


^  =  7777  E  (Wn>m)-m)2,  (38) 

A'fc  ( 

(n,m)G  Gfc 

and  |Cfc|  denotes  the  number  of  pixels  in  class  G),.. 

6)  To  avoid  sensitivity  to  resolution,  the  quantity  we  will  use  to  measure  the  improvements  to  the 
image  J  due  to  algorithm  cilg  is 

Qi.i,  =  ^  •  (39) 

We  refer  to  this  relative  performance  metric  Qialg  as  the  ultrasound  despeckling  assessment  index 
(USDSAI). 

If  the  classes  are  well  separated,  then  the  numerator  in  equation  (36)  will  be  large.  If  the  segmentation 
or  restoration  algorithm  produces  small  intra-class  variances,  then  the  denominator  of  equation  (36)  will 
be  small,  resulting  in  the  USDSAI  quantity  Qialg  to  be  large.  Large  value  would  indicate  that  alg  produces 
the  desirable  restoration  result.  If  Qialg  is  greater  than  one,  then  the  algorithm  being  tested  improved  the 
ultrasound  image  J.  Larger  Qialg  values  indicate  better  performance  when  comparing  various  algorithms. 
If  Qiaig  Is  less  If1311  one’  then  we  consider  the  algorithm  being  tested  as  detrimental  to  the  ultrasound 
image  J.  If  Qialg  is  equal  to  one,  then  it  is  suspect  that  any  improvements  were  accomplished. 

To  elucidate  the  significance  of  the  USDSAI  quantity  Qialg  defined  in  equation  (39),  we  examine  the 
three  extreme  cases.  The  first  case  is  the  simulated  image  J  is  not  constant  and  the  resulting  image  of 
some  algorithm  I aig  is  constant.  Clearly,  Qialg  =  which  we  define  as  zero.  Thus,  the  quantity  Qialg  is 
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equal  to  zero.  We  consider  algorithms  that  yield  small  Qialg  values,  less  than  one,  as  indicative  of  poor 
performance. 

In  the  second  case,  the  simulated  image  J  is  neither  a  constant  nor  a  scaled  version  of  the  template 
image.  The  resulting  image  of  some  algorithm  is  the  template  image,  that  is  I aig  =  T.  In  this  case 

K 

Qialg  =  J  =  00^  =  00  (40) 

when  the  average  of  the  squared  differences  of  the  inter-class  means  represented  as  K  is  not  equal  to  zero. 
We  will  consider  this  quantity  as  indicative  of  an  ideal  improvement  to  the  image  J  and  the  algorithm 
being  tested  performed  ideally. 

The  third  and  final  case  is  when  the  simulated  image  J  be  composed  of  classes  C\ ,  C'2,  C'3, . . . ,  CN 
and  the  pixel  values  of  the  resulting  image  I (dg  is  some  constant  except  for  a  single  point  within  each 
class  C\,  C'2,  C3, . . . ,  Cat,  without  loss  of  generality,  say  the  pixel  values  are  all  zero  except  at  differing 
points  in  each  class  have  values  Ki\Ri\,  iC2|i?2|,  K3\R3\, . . . ,  KN\RN\,  where  K,  >  0.  In  other  words,  the 
algorithm  being  tested  sets  the  pixels  in  each  class  to  the  same  constant  except  one  pixel  differs  from  the 
rest.  In  this  case,  the  sum  of  the  inter-class  mean  differences  squared  is 

N 

£(m  -  Vcf  =  £)(**  -  K,f  <  (N  -  1)  J2  Kl  (41) 

k>l  k>l  k= 1 

The  inequality  of  equation  (41)  is  attain  from  the  safe  assumption  that  >  0  for  all  k.  The  variance 
for  each  class  Ck  where  k  =  1,  2,  3, . . . ,  iV  is 

<4  =  (|C*|  -  1)K|.  (42) 

The  sum  of  the  intra-class  variances  is 

N  N 

X>a  =  £( \ct\-ml  (43) 

k=  1  k=l 

Let  Cmin  =  min{|Cfc|  where  k  —  1,  2,  3, ... ,  N},  so  that 

N  N 

k= 1  k= 1 


(44) 
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Using  the  inequalities  from  equations  (41)  and  (44),  the  quantity  Qialg  defined  in  equation  (36)  is  bounded 
above  by 

N 

(JV- 1)£>2 

Qi.„  < - r-  (45) 

_  U  fy)  if)  -L 

(iV -1)E*1 

fe=l 

In  this  present  case,  the  algorithm  being  evaluated  produces  an  image  composed  of  a  scaled  Kronecker 
delta  function  in  each  class.  When  Cmin  is  very  large,  the  original  unprocessed  image  contain  large  classes 
of  homogeneous  regions  and  the  metric  Qialg  returns  a  small  number.  Thusly,  the  Qialg  value  would  be 
smaller  than  Qj. 

To  present  an  objective  comparison,  we  perform  an  exhaustive  search  varying  the  parameters  of  each 
algorithm  so  that  the  USDSAI  value  Q[ ™  was  maximized.  We  used  these  optimal  parameters  of  each 
algorithm  to  perform  the  comparison  with  the  other  test  images.  The  result  of  this  exhaustive  search  for 
the  Wiener  filter  PSF  that  was  derived  from  the  simulation  is  a  two  dimensional  Gaussian  function  with 
horizontal  and  vertical  means  of  zero  and  standard  deviations  of  0.1. 

The  resulting  images  using  the  Nagao,  Lee,  Frost,  Kuan,  AWMF,  SRAD,  Wiener,  and  SBF  are  shown 
in  Fig.  2  and  3.  The  USDSAI  values  Q r  of  each  filter  is  given  in  Table  I.  The  observant  readers  will 

alg 

notice  that  the  subjective  results  shown  in  Fig.  2  and  3  reflect  the  quantitative  results  of  Table  I.  The 
average  of  all  USDSAI  values  Q\alg  for  all  algorithms  tested  is  given  in  the  rightmost  column  of  Table  I. 

The  resulting  images  produce  by  the  Nagao  filter  are  shown  in  Figs.  2(a)-  2(d).  A  close  inspection  of 
these  images  show  that  the  Nagao  filter  produces  small  “patchy”  regions  in  place  of  the  speckle.  More 
detrimental,  it  is  very  evident  that  the  edges  of  the  various  circular  regions  are  greatly  distorted.  The 
USDSAI  values  Of  and  Qpi  given  in  Table  I  are  less  than  one,  indicating  that  the  Nagao  filter 
is  counter  productive  in  improving  the  contrast  enhancement  of  these  simulated  ultrasound  images.  An 
intraclass  reduction  in  the  variances  of  the  different  classes  of  the  simulated  images  shown  in  Fig.  2(c)-  2(d) 
produced  USDSAI  values  Qfu  and  Qpv  ,  to  be  greater  than  one.  Thus,  a  overall  modest  quantitative 
improvement  due  to  the  Nagao  filter  is  indicated  by  the  average  USDSAI  value  slightly  greater  that  one 
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and  given  in  the  rightmost  column  of  Table  I. 

The  images  shown  in  Fig.s  2(e)-  3(d)  are  due  to  the  application  of  the  Lee,  Frost,  Kuan  filters, 
and  AWMF,  resp.  Compared  with  the  simulated  images  shown  in  Fig.  1,  the  subjective  or  qualitative 
improvements  due  to  these  filters  are  at  best  trivial,  if  any  improvements  can  even  be  claimed.  The 
USDSAI  values  of  these  four  algorithm  are  given  in  Table  I.  The  rightmost  column  of  Table  I  shows 
the  average  USDSAI  from  the  four  simulated  images  of  each  algorithm.  A  qualitative  visual  inspection 
and  quantitative  assessment  using  the  USDSAI  values  of  the  Lee,  Frost,  Kaun  filters  and  AWMF  indicate 
lackluster  results  when  compared  with  the  image  quality  and  USDSAI  values  of  SRAD,  Wiener  filtering, 
and  SBF. 

A  qualitative  assessment  of  the  resulting  images  after  applying  SRAD  is  given  in  Fig.  3(e)-  3(h)  and 
visibly  shows  that  the  variances  in  the  different  classes  are  decreased.  The  decrease  variances  in  the 
different  classes  due  to  SRAD  is  indicated  by  significantly  larger  USDSAI  values  shown  in  Table  I.  The 
rightmost  column  in  the  SRAD  row  of  Table  I  gives  the  average  USDSAI  value  and  is  noticeably  greater 
than  the  average  USDSAI  values  attained  by  the  Nagao,  Lee,  Frost,  Kuan  filters,  and  AWMF. 

The  result  of  Wiener  filtering  shown  in  Fig.  3(i)-  3(1).  Except  for  the  results  due  to  SBF,  these  results 
are  perceptively  better  than  the  images  produced  by  the  other  algorithms.  The  different  classes  of  each 
image  are  well  defined  and  the  edges  are  noticeably  more  pronounced  than  in  the  unprocessed  simulated 
images  shown  in  Fig.  1(e)-  1(h)  and  the  results  of  the  other  algorithms.  The  enhancements  provided  by 
Wiener  filtering  results  in  an  increased  USDSAI  value,  even  over  SRAD.  Unfortunately,  it  could  also 
equally  be  argued  that  a  low  pass  version  of  the  speckle  persists  and  reduction  in  the  intraclass  variance 
could  be  improved. 

The  resulting  four  images  produced  by  the  SBF  despeckling  method  using  a  applied  to  the  four  Field  II 
simulated  images  are  shown  in  Figs.  3(m)-  3(p).  It  is  visually  obvious  that  the  edges  of  the  various  disks 
are  equivalently  as  well  preserved  when  compared  with  the  results  of  the  Wiener  filter.  Equally  important 
the  intraclass  variance  is  decreased  in  each  simulated  image.  These  noticeable  contrast  improvements 
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(e)  J1 


(f)  J11 


(g)  3UI 


(h)  3IV 


Fig.  1.  Two  and  three  class  templates  and  ultrasound  simulations  used  in  the  comparison. 


are  reflected  by  the  USDSAI  values  in  the  last  row  of  Table  I.  Expect  for  Q v  ,  which  is  slightly 
less  than  Qji  ,  the  other  quantitative  results  show  improved  contrast  enhancement  over  the  other 
algorithms,  even  over  the  Wiener  filter.  The  average  USDSAI  values  of  the  SBF  method  indicate  greater 
contrast  improvements  over  the  Wiener  filter  and  significantly  greater  contrast  improvements  over  the 
other  despeckling  algorithms. 


B.  Using  Actual  Phantom  Image 

The  SRAD  and  SBF  despeckling  methods  that  performed  well  on  simulated  images  were  evaluated 
using  an  actual  ultrasound  scan  of  a  phantom.  Since  the  Wiener  filter  requires  prior  knowledge  of  the  PSF 
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Qi1 1 

alg 

Qi1] 

alg 

Qiin 

alg 

Qiiv 

alg 

Average 

F 

Nagao 

0.9901 

0.9730 

1.0315 

1.3651 

1.09 

I 

Lee 

1.2014 

1.1760 

1.2007 

1.2097 

1.197 

L 

Frost 

1.0479 

1.0211 

1.0446 

1.0457 

1.04 

T 

Kuan 

1.0007 

1.0006 

1.0012 

1.0016 

1.001 

E 

AWMF 

1.1005 

1.0936 

1.0959 

1.1025 

1.098 

R 

SRAD 

1.5344 

1.4729 

1.8603 

1.8187 

1.6716 

S 

Wiener 

2.0884 

1.8423 

2.0256 

2.0072 

1.9909 

SBF 

1.9974 

2.1160 

2.1491 

2.1144 

2.0942 

TABLE  I 

USDSAI  VALUES  Qji  FOR  THE  VARIOUS  ALGORITHMS  TESTED. 


and  this  paper  is  not  focused  on  estimating  or  calibrating  the  PSF  of  the  ultrasound  system,  the  Wiener 
filter  was  not  evaluated.  Rather,  we  visually  evaluate  the  performance  of  SRAD  and  SBF  at  various 
iterations.  This  evaluation  uses  an  actual  B  mode  image  taken  from  a  scan  of  a  real  phantom.  The  original 
phantom  image  is  shown  in  the  top  image  of  Fig.  4.  The  phantom  consist  of  three  equal  size  disks  with 
varying  brightness.  The  three  images  in  Figs.  4(b),  4(d),  and  4(f)  show  the  results  of  SRAD  at  10,  25, 
and  50  iterations,  resp.  The  images  on  the  right  column  of  Fig.  4  show  in  Figs.  4(c),  4(e),  and  4(g)  are 
the  results  of  the  SBF  method  using  an  11  x  11  window  at  10,  25,  and  50  iterations,  resp.  All  the  images 
in  Fig.  4  are  shown  scaled  to  eight  bit  dynamic  range,  that  is  256  gray  levels.  A  visual  comparison  of 
these  images  show  that  at  the  three  iterations  the  brightness  of  rightmost  disk  is  greater  than  the  middle 
disk  and  the  brightness  of  the  middle  disk  is  greater  than  the  the  leftmost  disk  with  the  results  of  the 
SBF  method  than  SRAD.  Additionally,  the  leftmost  darkest  disk  at  50  iterations  almost  blends  into  the 
background,  while  even  at  the  50th  iteration  of  SBF  this  disk  is  well  pronounced.  It  is  well  worth  noting 
that  the  bright  dots  between  the  leftmost  and  the  middle  disks  are  better  preserved  with  SRAD  than  SBF. 


From  the  result  of  this  evaluation  using  an  actual  ultrasound  phantom  image,  we  conclude  that  SBF  is 
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better  at  preserving  the  contrast  of  the  three  disks  at  a  cost  of  decimating  the  fine  resolution  features  such 
as  the  bright  dots  between  the  leftmost  and  middle  disks. 

C.  3-D  GVF  Active  Contour  Surface  Rendering 

Preserving  or  enhancing  the  contrast  of  large  blunt  objects  at  the  cost  of  losing  or  fading  smaller  objects 
is  preferred  in  our  ultimate  goal  of  rendering  a  3-D  surface  from  2-D  ultrasound  slices  of  an  organ  such  as 
a  heart,  prostate,  kidney,  etc.  In  this  evaluation  we  used  an  linear  array  transducer  to  scan  an  egg  shaped 
phantom.  The  2-D  slices  are  taken  at  every  1mm  along  the  long  axis  of  the  egg  phantom.  There  is  a  total 
of  41  slices  and  the  B  mode  middle  slice  of  the  egg  phantom  is  shown  in  Fig.  5(a).  The  3-D  rendering 
using  a  naive  implementation  of  SRAD  processed  slices  yielded  a  surface  that  in  no  way  resemble  the  egg 
phantom.  So  each  slice  was  SRAD  processed  by  meticulously  adjusting  the  parameters  and  iterations  of 
SRAD  so  that  a  well  defined  high  contrast  image  was  produces.  After  each  slice  was  optimally  visually 
contrast  enhanced  by  SRAD,  a  post  despeckling  algorithm  of  morphologically  filling  in  the  holes  was 
applied.  The  result  of  the  optimized  SRAD  followed  by  morphologically  filling  the  holes  applied  to  the 
image  in  Fig.  5(a).  Albeit  labor  intensive,  the  result  of  these  processes  are  exceptionally  excellent.  We 
will  used  the  3-D  GVF  active  contour  rendering  of  these  exceptional  SRAD  with  morphological  hole 
filling  slices  as  ground  truth,  a  gold  standard  to  compare  the  results  of  processing  with  SBF.  We  applied 
the  same  75  iteration  SBF  process  using  a  11  x  11  window  to  each  of  the  41  slices.  This  processing  was 
implemented  using  Mathworks  Matlab  7.0.0.19920  on  a  Intel  P4  3.8  GHZ  CPU.  The  processing  time 
need  by  the  SBF  to  process  all  41  slices  was  8.8268  minutes.  The  result  of  SBF  applied  to  the  middle 
slice  is  shown  in  Fig.  5(c). 

The  3-D  rendering  of  the  final  contour  using  the  same  number  of  iteration  of  the  3-D  GVF  snake 
produced  from  the  41  slices  processed  by  the  human  optimized  SRAD-morphological  hole  filling  operation 
and  SBF  are  shown  in  Figs.  6(a)  and  6(b),  resp.  The  3-D  renderings  are  displayed  so  that  the  worst 
part  of  each  3-D  surface  is  shown.  In  the  3-D  rendering  provided  by  the  human  optimized  SRAD  with 
morphological  hole  filling,  there  is  an  erroneous  indention  just  above  the  horizontal  bisecting  plane.  The 
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3-D  rendering  using  the  SBF  processed  slices  exhibits  two  unwanted  steps  above  the  horizontal  bisecting 
plane.  Both  3-D  GVF  active  contour  rendering  using  the  human  optimized  SRAD  with  morphological 
holes  filled  and  SBF  processed  perform  well  in  that  overall  shape  of  the  egg  is  captured.  Estimation  of  the 
volume  enclosed  by  each  the  surface  created  from  the  human  optimized  SRAD  with  morphological  hole 
filling  errs  by  approximately  7%  of  the  true  volume.  Volume  estimate  enclosed  by  the  3-D  GVF  surface 
using  SBF  processed  slices  yields  a  10%  error  to  the  true  volume.  The  slightly  high  error  using  the  3-D 
rendering  from  the  SBF  processed  slices  is  insignificant  when  one  considers  the  effortlessness  of  attaining 
SBF  processed  slices  over  optimizing  each  slice  with  SRAD  and  subsequently  applying  a  morphological 
hole  filling  operation. 

V.  Conclusion 

The  removal  or  reduction  of  speckle  noise  while  preserving  or  enhancing  edge  information  of  an 
ultrasound  image  is  an  extremely  difficult  task.  Since  the  GVF  is  derived  from  an  edge  map,  robustly 
despeckling  each  slice  is  vital  prior  to  apply  the  3-D  GVF  snake  to  determine  a  surface  from  ultrasound 
slices.  We  consider  of  a  wide  variety  of  filtering  algorithms  for  this  pre-rendering  step.  The  despeckling 
methods  describe  in  this  paper  and  in  general  scrutinize  every  pixel  values.  We  present  a  novel  iterative 
despeckling  method  SBF  that  at  each  iteration  only  scrutinize  the  value,  which  are  considered  outliers. 
Without  using  the  outlying  values  to  determine  the  local,  the  SBF  method  replaces  these  outliers  by  values 
that  approach  the  mean  of  the  local  homogeneous  region.  Each  iteration  of  this  new  method  compress 
the  image  pixel  values  so  that  the  differences  in  interclass  means  are  preserved  or  possibly  enhanced 
while  the  intraclass  variance  is  decreased.  The  superior  contrast  enhancement  preformance  of  the  SBF 
method  is  established  by  our  experimentation  using  Field  II  simulated  ultrasound  images  and  evaluated 
with  USDSAI  performance  metric  defined  in  equation  (39).  Our  next  experimentation  evaluation  compare 
the  anisotropic  diffusion  method  SRAD  with  SBF  on  an  actual  ultrasound  image  of  a  phantom  consisting 
of  three  disks  of  the  same  size  but  varying  brightness.  The  SBF  exhibited  excellect  contrast  improvements 
of  the  three  disks  over  SRAD,  but  at  the  cost  of  fading  the  small  bright  spots.  This  tradeoff  is  preferred 
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in  our  ultimate  goal  of  applying  a  3-D  GVF  snake  to  render  an  accurate  surface.  Our  last  experiment 
shows  that  the  surface  found  by  the  3-D  GVF  snake  using  the  SBF  processed  slices  is  comparable  in 
size,  shape,  and  volume  of  the  surface  rendered  by  applying  the  3-D  GVF  snake  to  the  tediously  human 
optimized  SRAD  and  morphologically  hole  filled  slices. 
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Fig.  2.  The  resulting  images  from  various  adaptive  filtering  algorithms. 


II 

IRAD 


W  iener 


Fig.  3.  The  resulting  images  from  AWMF.  SRAD.  Wiener,  and  SBF  filters. 
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(d)  25  iterations  of  SRAD. 
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(c)  10  iteration  of  SBF. 
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(e)  25  iteration  of  SBF. 

• 

(f)  50  iterations  of  SRAD. 


(g)  50  iteration  of  SBF. 


Fig.  4.  The  resulting  images  from  10,  25,  50  iterations  of  SRAD  and  SBF. 
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followed  by  a  morphological  hole 
filling  processed  middle  slice. 


(c)  The  SBF  processed  middle 
slice. 


Fig.  5.  The  original  middle  slice,  the  processed  slice  using  SRAD  (optimized  by  human  interactions)  and  a  morphological  hole  filling 
operation,  and  SBF  using  an  11  x  11  window. 
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(a)  3-D  surface  rendering  of  an  egg  using  a  human  (b)  3-D  surface  rendering  of  an  egg  using  SBF. 

optimized  SRAD  followed  a  morphological  hole  filling 

operation. 

Fig.  6.  The  3-D  surface  renderings  of  the  egg  phantom  where  each  slice  is  processed  with  the  visually  optimized  SRAD  followed  by  a 
morphological  hole  filling  operation  and  the  SBF  using  an  11  x  11  window. 
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(c)  The  top  view  of  egg  using  SRAD. 


(d)  The  top  view  of  egg  using  SBF. 


Fig.  7. 


Side  and  top  views  of  the  3-D  rendering  of  the  egg  phantom  using  the  human  optimized  SRAD  followed  by  a  morphological  hole 


filling  operation  and  the  surface  rendering  using  SBF. 
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Abstract 

This  paper  proposes  a  robust  method  to  restore  piecewise  constant  signals  contaminated  by  additive 
and/or  multiplicative  noise.  Additionally,  the  signal  maybe  distorted  by  the  underlying  physics  of  the 
signal  acquistion  system.  The  possible  distortion  caused  by  an  signal  acquistion  system  is  modeled  as  a 
convolution  with  some  impulse  response.  Analytic  analysis  will  establish  that  when  the  signal  is  view 
as  a  stochastic  process  and  regardless  of  whether  the  noise  is  additive  or  multiplicative,  exact  or  perfect 
restoration  is  pausible.  The  implementation  of  these  well  known  and  firmly  established  analytic  facts 
about  random  processes  will  utilize  a  stochastically  driven  filter  that  we  appropriately  call  the  squeeze 
box  filter.  We  will  empirically  show  that  the  squeeze  box  filter  can  provide  a  faithful  and  accurate 
solution  to  the  additive  and  multiplicative  noise  restoration  problems  under  the  regularization  constraint 
that  the  ideal  or  desired  signal  is  piecewise  constant  and  provided  that  the  expected  value,  i.e.  mean 
or  average,  of  the  stationary  noise  is  known  a  priori.  Even  if  the  expected  value  of  the  noise  is  not 
known,  the  proposed  filtering  method  produces  a  robust  method  to  decrease  variances  within  intervals 
where  the  signal  is  a  constant  value,  while  the  opposing  task  of  preserving  edges  is  maintained. 
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A  Stochastic  Approach  to  the  Noise 
Restoration  Problems 

I.  Introduction 

The  physical  nature  of  ultrasonic  sound  wave  as  an  imaging  modality  allows  a  safe  (uses 
non-ionizing  energy),  noninvasive,  and  relatively  inexpensive  medical  diagnostic  imaging  tool. 
Though  these  justifications  for  using  ultrasound  images  as  a  medical  diagnostic  tool  are  ap¬ 
pealing,  the  images  provided  by  an  ultrasonic  imaging  system  is  relatively  unfavorable  when 
compared  to  systems  that  use  X-rays  (e.g.  film-screen  or  digital  radiography,  fluoroscopy,  com¬ 
puted  tomography)  and  nuclear  magnetic  resonance.  The  basic  physics  of  ultrasound  imaging 
and  other  medical  imaging  modalities  are  well  described  in  [1], 

As  with  all  known  nonbiological  imaging  systems,  the  images  produced  by  an  ultrasonic 
imaging  system  are  degraded  by  some  form  of  noise.  The  noise  inherent  to  ultrasound  imaging  is 
dependant  on  several  factors  such  as  axial/lateral  resolutions  of  the  transmit/receive  frequencies, 
axial  and  lateral  distances  from  the  transducer,  the  image  resolution,  the  characteristics  of 
the  tissue(s),  organ(s),  or  object(s)  being  scanned,  electronic  noise,  scalar  quantization  of  the 
reflectivity  values,  etc.  Though  the  appearance  of  the  noise  can  vary,  the  inherent  noise  can 
corrupt  the  ultrasound  system  by  rendering  an  image  with  a  grainy  texture,  that  is  the  image 
appear  contaminated  with  “salt  and  pepper”  like  texture.  This  type  of  noise  is  commonly  called 
speckle  by  the  ultrasound  and  synthetic  aperture  radar  communities.  To  show  how  speckle 
appears  in  a  typical  medical  ultrasound  image,  we  offer  Fig.  1(a),  which  displays  a  B-mode 
logarithmically  compressed  envelop  detected  ultrasound  image  of  a  phantom.  The  phantom 
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consist  of  three  disks  each  with  different  intensity  values  and  are  approximately  50  pixels  in 
diameter.  The  centers  of  each  disk  are  approximately  on  row  225.  The  phantom  also  contains 
several  point  scatters  with  large  intensities  along  row  90  and  along  column  120.  The  noise  is 
evident  in  this  image  where  the  graininess  is  prominent  throughout  the  entire  image  and  the 
salt  and  pepper  texture  could  appear  coarser  in  the  far  field1  (not  evident  in  this  image).  Also 
evident  in  Fig.  1(a)  are  the  shadow  artifacts  that  occur  away  from  the  transducer  immediately 
below  each  high  intensity  point  scatters  and  appear  like  darkened  comet  tails.  Though  not  as 
visually  evident,  a  shadow  artifact  is  present  below  the  center  disk.  This  shadowing  phenomena 
is  a  known  problem  with  ultrasound  imaging  and  results  in  smaller  intensity  values  on  average 
within  the  shadow  regions  to  erroneously  occur.  The  shadowing  problem  is  beyond  the  scope 
this  paper  and  we  will  address  this  in  our  future  research. 

In  Fig.  1(b)  we  show  the  intensity  profile  of  row  225  of  Fig.  1(a).  This  row  approximately 
bisects  all  three  disks  throught  the  center.  From  the  profile  in  Fig.  1(b)  it  is  difficult  to  determine 
which  intervals  corresponds  to  the  various  three  disks.  Even  more  difficult  is  to  determine  exactly 
the  actual  reflectivity  values2  of  each  phantom  disk,  since  the  processing  used  to  produce  a 
visually  appeasing  image  skews  these  values.  In  Fig.  1(c)  is  a  plot  of  a  possible  noise  free 
signal  of  the  intensity  profile  signal  shown  in  Fig.  1(b).  This  piecewise  constant  signal  was 
determine  by  taking  the  average  in  intervals  that  are  known  to  be  a  constant.  It  can  be  observed 
from  the  proposed  noise  free  signal  shown  in  Fig.  1(c)  that  the  intensity  value  in  interval 
[50, 100]  corresponds  to  the  left  most  black  (low  intensity  value)  disk,  the  intensity  value  in 
interval  [140, 190]  corresponds  to  the  middle  disk,  and  the  intensity  value  in  interval  [240,  290] 

'The  region  away  from  the  transducer  and  pass  the  focal  point. 

“The  actual  reflectivity  value  directly  relates  to  tissue/organ  elasticity,  compressibility,  and/or  density. 
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corresponds  to  the  right  most  brightest  (high  intensity  value)  disk  of  Fig.  1(a). 

In  Fig.  1(d),  we  show  the  intensity  profile  of  column  165,  which  approximately  vertically 
bisects  the  middle  disk  of  Fig.  1(a).  Fig.  1(e)  is  a  plot  of  a  piecewise  constant  signal  determined 
by  taking  the  average  with  intervals  that  are  known  to  be  constant.  The  lower  intensity  value 
in  the  right  most  interval  of  the  plot  in  Fig.  1(e)  is  a  shadowing  artifact.  Ideally,  the  average 
intensity  value  in  this  interval  should  be  the  same  as  the  average  intensity  value  in  the  left  most 
constant  interval  of  Fig.  1(e). 

In  creating  the  noise  free  piecewise  constant  signal  in  Fig.s  1(c)  and  1(e),  we  used  our 
knowledge  of  the  phantom  to  determine  the  values  and  location  of  each  constant  intervals, 
that  is  we  knew  where  the  edges  are  located.  This  in  all  practical  and  clinical  use  of  medical 
ultrasound  imaging  as  a  diagnostic  tool  is  not  the  case.  Determining  a  possible  noise  free  signal 
such  as  the  ones  shown  in  Fig.s  1(c)  and  1(e)  from  the  noisy  signals  shown  in  Fig.s  1(b)  and  1(d) 
is  a  classical  signal  restoration  problem. 

The  restoration  problem  of  determining  the  ideal  signal  from  a  distorted  and  noisy  signal  is 
considered  ill  posed  in  which  the  solutions  are  nonunique  and/or  it  maybe  difficult  to  verify 
the  significance  of  the  proposed  solution.  To  consider  the  problem  well  posed,  either  a  priori 
knowledge  of  the  desired  signal  is  needed  or  constraints  on  the  ideal  signals  are  required.  The 
latter  is  called  regularization.  Some  regularization  constraints  for  signals  are  proposed  in  [2],  [3]. 
Some  examples  of  two  dimensional  regularization  constraints  for  images  are  proposed  in  [4]- 
[7].  To  avoid  being  clouded  by  dimensionality,  the  explanations  and  examples  given  in  this 
paper  are  strictly  using  one-dimensional  signals.  This  allows  us  to  present  our  ideas  on  signal 
restoration  in  a  clear  concise  manner  without  being  bothered  by  the  technicalities  of  extending 
one-dimensional  algorithms  to  multiple  dimensions.  Since  the  composition  of  ultrasound  images 
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is  an  synchronous  arrangement  of  one  dimensional  signals,  the  analysis  and  processing  of  one 
dimensional  signal  in  this  paper  is  relevant  to  enhancing  and  improving  ultrasonic  imaging. 

The  mathematical  modeling  of  ultrasound  data  is  debateable.  The  evidence  provided  by  Wagner 
et.  al.  in  [8],  [9]  suggest  that  ultrasound  speckle  is  multiplicative  and  Rayleigh  distributed.  While 
Loupas  et.  al.  in  [10]  and  Karaman  et.  al.  in  [11]  claim  speckle  is  composed  of  a  combination 
of  multiplicative  and  additive  components3.  The  determination  of  the  best  method  to  model 
ultrasound  data  is  beyond  the  scope  of  this  paper.  Instead  this  paper  will  focus  on  the  usual 
additive  and  multiplicative  signal  restoration  problems. 

The  two  usual  restoration  problems  [12]  are  to  determine  the  ideal  signal  f(n)  from  the 
detected  signal  g(n)  given  by  the  additive  model 

g(n)  =  p(n)  *  f(n )  +  77(71)  (1) 

or  by  the  multiplicative  model 

g(n)  =  p(n)  *  {f(n}g(n)}  (2) 

where  p(n)  denotes  the  impulse  response  of  some  system,  the  asterick  symbol  represents  linear 
convolution,  and  77(71)  is  a  random  variable  with  some  probability  density  function  (PDF).  The 
noise  samples  77(77)  are  generally  assumed  to  be  white  and  Gaussian  distributed.  Moreover, 
the  noise  is  typically  assumed  to  be  independent  of  f(n).  Generally,  there  is  an  unmentioned 
assumption  that  the  noise  is  a  wide  sense  stationary  process  and  we  will  accept  this  underlying 
assumption  throughout  this  paper.  In  the  models  given  in  equations  (1)  and  (2),  the  distorted 
signal  is  modeled  as  a  convolution  by  an  impulse  response,  79(77).  This  distortion  is  usually  a 

’They  do  not  mention  the  noise  distribution. 
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blurring  or  smoothing  of  the  desired  signal.  Thus  the  impulse  response  is  in  nature  a  low  pass 
filter. 

The  additive  restoration  problem  given  in  equation  (1)  has  extensively  been  studied  and  many 
solutions  to  this  problem  have  been  employed  with  varying  success.  These  solutions  include 
linear  filtering  such  as  a  moving  average,  non-linear  filters  such  as  the  various  adaptive  median 
and  mean  filters  [4],  [23],  Weiner  filtering  [13],  wavelet  thresholding  [12],  etc. 

Methods  to  solve  the  multiplicative  restoration  problem  of  equation  (2)  are  and  have  been  of 
great  interest  to  ultrasound  and  SAR  image  processing  societies.  Solutions  to  the  two  dimensional 
extension  of  the  multiplicative  noise  model  have  been  proposed  in  [7],  [10],  [11],  [14]-[22].  None 
of  the  proposed  method  especially  for  multiplicative  noise  restoration  claim  an  exact  restoration 
of  the  ideal  signal  and  may  only  offer  subjective  improvements.  Provided  the  usual  assumptions, 
the  noise  is  stationary  and  is  a  zero  mean  Gaussian  distributed  random  process,  are  accepted, 
exact  reconstruction  of  the  ideal  signal  will  be  shown  to  be  analytically  possible  from  both 
the  additive  model  given  in  equation  (1)  or  the  multiplicative  model  given  in  equation  (2). 
Though  exact  reconstruction  is  possible  in  an  analytic  sense,  the  practical  implementation  of 
these  analytic  concepts  is  not  without  trade-offs.  Namely,  the  task  of  determining  the  charac¬ 
teristics  of  a  random  process,  e.g.  mean  and  variance,  from  only  a  few  recorded  samples  is 
debateable.  Nevertheless,  we  will  reasonably  compromise  when  required  and  develop  a  practical 
implementation  of  the  analytic  concepts  presented  in  this  paper.  Our  examples  will  empirically 
show  using  the  assumptions  that  the  noise  is  white,  normally  distributed,  and  stationary  that  the 
stochastically  driven  proposed  filter  robustly  extracts  the  ideal  piecewise  constant  signal. 
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II.  Various  Signal  Restoration  Methods 

In  this  section  we  will  describe  a  variety  of  signal  restoration  methods.  The  performance  of  our 
proposed  method  will  be  compared  with  the  performance  of  the  filters  described  in  this  section. 
The  quantitative  and  qualitative  assessment  of  this  comparison  will  be  given  in  section  V. 

A.  Order  Statistic  Filters 

In  [23]  Bovik  et.  al.  described  a  family  of  order  statistic  filters  (OSF).  The  OSF  replaces  each 
and  every  sample  of  the  noisy  signal  with  the  weighted  sum  of  ordered  values  taken  from  an 
odd  length  window.  More  precisely,  given  a  window  of  odd  length  N,  the  value  of  the  current 
sample  is  set  to 

n+M 

9osf{ti )  =  ^2  (3) 

i=n—M 

where  M  =  [yj,  [•]  is  the  greatest  integer  function,  7j(i)  e  {g(i)  \  n  —  M  <  i  <  n  +  M},  and 
g(i)  are  ordered  so  that  g(i)  <  g(i  +  1)  for  all  i  <G  [n  —  M,n  +  M  —  1],  These  weights  cu 
are  computed  so  that  the  expected  value  of  the  filtered  signal  is  optimal  in  the  mean  squared 
error  sense.  The  optimally  of  the  OSF  is  under  the  constraint  that  the  signal  is  constant  and  the 
weights  of  the  filter  are  dependant  on  the  additive  noise  distribution.  For  example,  if  oi%  =  -7  for 
all  i  e  [n  —  M,n  +  M],  then  the  filter  described  in  equation  (3)  is  the  moving  average  filter4. 
Bovik  et.  al.  analytically  show  the  moving  average  is  optimal,  that  is  the  expected  value  of  the 
mean  squared  error  is  minimized,  when  the  ideal  signal  is  constrainted  as  an  constant  signal  and 
the  noise  is  zero  mean  and  normally  distributed.  If  an  =  1  and  cu  =  0  for  all  i  ^  n,  then  the 
filter  in  equation  3  is  the  median  filter.  The  justifications  in  [23]  does  not  claim  that  the  median 

4The  moving  average  filter  is  technically  not  a  OSF,  since  ordering  of  the  windowed  values  are  not  required.  We  place  this 
example  in  this  section  purely  for  the  purpose  of  efficiency. 
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is  optimal  with  any  of  the  considered  noise  distribution.  We  will  include  the  median  filter  in  our 
comparison  in  section  V  so  that  a  wide  range  of  OSF  filters  are  represented  in  the  experimental 
comparison.  The  last  OSF  we  will  include  in  our  comparative  study  has  a  window  length  of 
nine  and  the  following  symmetric  weight  values: 


Oin 

0.36469 

C^n—1 

0.23795 

C^n+1 

&n— 2 

0.06965 

—  ^n+2 

3 

0.02904 

^71+ 3 

&n— 4 

-0.01899 

C^n+4 

In  [23]  Bovik  et.  al.  claims  the  OSF  with  the  weights  given  in  equation  4  is  optimal  in  minimizing 
the  expected  value  of  the  mean  squared  error  for  a  Laplacian  noise  distribution  and  provided 
that  the  ideal  signal  is  a  constant. 


B.  Wiener  Filter 


The  Wiener  filter  [13]  is  a  frequency  domain  solution  to  the  additive  noise  problem  in 
equation  1.  Let  g(n)  < — »  G  (e?"’  )  be  discrete  time  Fourier  transform  (DTFT)  be  pairs,  then 
the  Wiener  filter  is  implemented  in  the  Fourier  domain  as 

gWnr(n)  <—  Gwnr  (e^)  =  G  (ej“)  W  (e?w)  (5) 


where 


W  (ejuJ) 


P*  ( ejuJ ) 

P  (e^)  |2  +  <x2 


(6) 


and  cr2  is  the  variance  of  the  additive  noise  component.  It  is  important  to  note  that  the  performance 


of  the  Wiener  filter  is  dependant  on  a  priori  knowledge  or  an  approximation  of  the  transfer 


function  of  p(n)  and  the  noise  variance  cr2. 
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C.  Wavelet  Thresholding 

It  has  been  shown  in  [24]  and  [25]  that  applying  a  threshold  operator  to  the  discrete  wavelet 
transform  (DWT)  coefficients  is  a  useful  method  to  remove  or  reduce  additive  white  noise.  The 
DWT  of  a  signal  described  by  Mallat  in  [26]  decomposes  the  signal  into  various  subbands. 
Assuming  that  the  noise  power  is  below  some  threshold,  then  applying  a  hard  [25]  or  a  soft 
thresholding  [24]  operation  to  the  DWT  coefficients  reduces  the  noise  power.  In  effect  the  inverse 
DWT  (IDWT)  reconstruction  of  the  threshold  coefficients  produces  an  approximate  noise  free 
signal. 

Let  w(n)  =  DWT{g(n)}  be  the  DWT  of  g(n).  For  T  >  0  the  hard  thresholding  operation  is 
defined  as 

I  w(n)  if  \w(n)\  >  T 

Whard{n )  =  <  (7) 

0  otherwise. 

The  soft  thresholding  operation  is  defined  as 

' 

w(n )  —  T  if  \w(n)\  >  T  and  w{n)  >  0 

wsoft(n )  =  <  w(n)  +  T  if  \w(n)\  >  T  and  w(n)  <  0  (8) 

w(n)  otherwise. 

An  approximation  of  the  additive  noise  free  ideal  signal  is  taken  to  be  the  IDWT  of  Whard(p)  or 
Wsoftiji).  The  robustness  of  the  thresholding  operations  are  dependant  on  the  choice  of  T  and 
a  optimal  choice  [24],  [25]  is  T  =  ori\j2  log  N  where  an  is  the  noise  standard  deviation  and  N 
is  the  length  of  g(n). 

III.  A  Stochastically  Driven  Filter 

Even  without  any  a  priori  knowledge  of  the  noise  characteristics  besides  stationarity,  the 
proposed  stochastically  driven  filter  will  provide  a  meaningful  restoration  of  a  piecewise  constant 
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signal.  We  use  the  term  “ meaningful ”  in  the  sense  that  the  proposed  restoration  method  will 
produce  a  constant  value  in  intervals  where  the  ideal  signal  is  a  constant  value.  Equally  as 
important,  the  opposing  task  of  preserving  step  up  and  step  down  edges  will  be  maintained. 

For  a  random  variable  X  it  is  well  known  and  can  be  shown  that  even  without  prior  knowledge 
of  the  stationary  noise  distribution  that 

E{cX  +  b}  =  cE{X}  +  b  (9) 

/OO 

xPrx(x)dx ,  and  Pi'  y(.x)  is  the  probability  density  function 

-OO 

(PDF)  of  the  random  variable  A".  Importantly,  determining  the  expected  value  or  mean  of  a 
random  value  will  preserves  edges,  that  is  for  constants  c i  ^  C2,  Xt  and  X2  are  independent 
and  identically  distributed  random  variables,  it  can  be  deduced  that 


E{clXl  -  c2x2}  =  (ci  -  c2)E{ X}  ^  0 


provided  that  E{X}  ^  0. 

From  the  well  know  fact  given  in  equation  (9)  and  without  loss  of  generality  let  b 
simple  algebraic  malipulation  yields 


_  E{cX} 
E{X}' 


0,  a 

(10) 


Provided  that  the  expected  value  of  the  random  variable  X  is  not  equal  to  zero,  the  triviality  of 
equation  (10)  yields  the  a  priori  unknown  value  c  provided  that  E{cX}  and  E{ X}  are  known 
or  can  be  determined.  In  the  context  of  the  multiplicative  noise  restoration  problem  modeled  in 
equation  (2)  and  ignoring  the  impulse  response5  p(n)  for  now,  we  arrive  at  the  nth  sample  of 
the  ideal  signal  from 


fH 


E{g(n)} 

E{v(n)} 


(ID 


5Restrict  the  impulse  response  to  be  the  delta  function.  We  will  counter  effect  the  impulse  response  later  in  this  paper. 
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where  g(n)  and  77(77)  are  cast  as  random  variables.  If  the  expected  value  of  the  noise  term  77(71) 
is  known  and  not  equal  to  zero,  then  to  achieve  exact  restoration  it  is  only  required  to  determine 
the  expected  value  of  the  detected  or  acquired  signal  #(77). 

To  determine  the  expected  value  of  g(n),  we  view  the  local  extrema,  that  is  the  local  peaks 
and  valleys,  to  be  outliers  of  some  PDF.  In  opposition  to  the  methods  in  section  II,  where  every 
sample  is  scrutinized,  the  filtering  method  presented  in  this  section  only  scrutinizes  the  samples 
that  are  considered  outliers.  These  outliers  are  replace  with  a  value  that  occurs  with  greater 
probability.  To  ensure  that  the  outliers  are  replaced  with  a  value  with  a  greater  probability,  we 
replace  them  with  the  mean  determined  from  some  local  neighborhood  J\f.  It  is  important  to  keep 

in  mind  that  the  local  extrema  are  considered  outliers  and  should  not  be  used  in  the  determination 

of  the  local  mean.  The  choice  of  the  neighborhood  M  is  extremely  important,  since  the  mean  of 
some  PDF  is  determined  by  samples  in  J\f.  This  produces  another  sequence  with  locally  reduced 
variance.  This  new  sequence  may  contain  local  extrema,  which  we  consider  as  outliers.  Thus,  the 
filtering  process  is  iterated  until  the  limiting  sequence  in  the  Cauchy  sense,  i.e.  the  root  signal, 
is  attained.  More  precisely,  to  arrive  at  E{g(n)}  for  all  77,  we  propose  the  following  iterative 
method. 

1)  Each  iteration  i  begins  by  determining  the  set  of  locations  of  local  maxima  (peaks)  and 
local  minima  (valleys).  The  locations  of  these  extrema  are  defined  by  the  set 

Me  =  {n  |  #7-1(77)  meets  condition  1  or  2  }  (12) 

Condition  1:  #7-1(77)  >  <77-1(77  —  1)  and  <77-1(77)  >  <77-1(77  +  1), 

Condition  2:  #7-1(77)  <  #7-1(77  -  1)  and  #7-1(77)  <  #7-1(77  +  1). 

The  local  peaks  and  valleys  of  a  length  twenty  randomly  generated  sequence  are  shown 
in  Fig.  2  as  A  and  y,  resp.  This  length  twenty  sequence  only  has  seven  samples  that  are 
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not  considered  a  local  extremum.  They  occur  at  sample  0,  3,  4,  6,  14,  18,  and  19.  The  first 
and  last  samples  are  never  considered  as  extrema  and  get  replaced  by  the  local  mean  after 
convergence  is  verified.  All  the  other  samples  are  either  a  local  peak  or  a  local  valley,  so 
in  this  example  the  set  Me  =  {1,  2, 5,  7,  8,  9, 10, 11, 12, 13, 15, 16, 17}. 


Fig.  2.  Step  1  of  SBF:  A  indicates  a  local  maximum  and  v  indicates  a  local  minimum. 


2)  Without  using  the  local  extrema  values,  poll  neighboring  samples  to  determine  the  local 
mean.  These  extrema  are  replaced  with  the  local  mean  values  that  is  for  n  G  Me 


9i(n) 


1 


9i- iM 

m£j\f 


(13) 


where  M  is  some  local  neighborhood  of  n,  |A/”|  is  the  cardinality  of  set  M,  and  n  (f  M. 
In  Fig.  3,  we  illustrate  with  astricks  (*)  the  possible  values  used  to  calculate  the  local 
means  of  a  local  maximum  (A)  and  a  local  minimum  (y)  in  Fig.s  3(a)  and  3(b),  resp. 
The  values  shown  as  A  and  y  are  not  used  to  calculate  the  local  mean.  In  Fig.  4  we  show 
the  replacement  values  of  each  extrema  as  “o”s. 
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(a)  A  local  maximum  is  shown  with  a  A. 
The  astricks  (*)  show  possible  values  used 
to  compute  the  local  mean. 


(b)  A  local  minimum  is  shown  with  a  V- 
The  astricks  (*)  show  possible  values  used 
to  compute  the  local  mean. 


Fig.  3.  Step  2.  The  *  denote  the  possible  sample  values  used  to  determine  the  local  mean.  Note  that  the  local  peak  shown  as 
a  A  and  the  valley  shown  as  v  are  not  used  to  determine  the  local  mean. 


3)  If  convergence  in  the  Cauchy  sense  is  not  attained,  that  is 

^2  \di-i(n)  ~  9i(n) |  >  e  (14) 

n 

for  some  predefined  e  >  0,  then  another  iteration  is  performed.  If  Cauchy  convergence 
is  attained,  then  g^n)  is  termed  the  root  signal.  Essentially  nonconvergence  indicates  a 
substantial  amount  of  peaks  and  valleys  still  exists.  Convergence  indicates  that  the  outlier 
values  are  removed.  The  ideal  signal  is  the  root  signal  in  the  additive  noise  case 

f(n)  ~  Qi{n)-  (15) 

In  the  multiplicative  noise  restoration  case  the  ideal  noise  free  signal  is  taken  to  be  the 
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Fig.  4.  Step  3  of  SBF:  local  mamimum  (A)  and  local  minimum  (v)  are  replaced  by  the  mean  within  a  local  neighborhood 
(°). 


root  signal  scaled  by  the  recipocal  of  the  nonzero  expected  value  of  the  stationary  noise, 


f(n) 


9i\n) 


E{v(n)}' 


(16) 


To  show  that  this  method  produces  a  convergent  sequence  of  sequences,  let  gi(n)  be  the 
sequence  produced  by  the  ith  iteration  of  the  proposed  method,  where  g0(n )  =  g(n).  Denote  the 
local  minimum  value  in  a  neighborhood  Af  of  sequence  gdn)  as  nrii  =  min {//,(//  )  j  nfAf}  and 
analogously  the  local  maximum  value  in  Af  of  the  sequence  fj,(n)  as  Mt  =  max{y,:(n)  |  n  G 
Af}.  It  is  clear  by  the  definition  that  rri,  <  gi(n)  <  Mi  or  equivalently  restated  g,  (n)  G  [m, ,  M{\ 
for  all  n  G  Af .  Suppose  at  some  iteration  i  that  g,(nr)  is  a  local  extremum,  we  replace  //,(///) 
with  the  mean  of  {giin)  \  n  G  Af  and  n  f  n'}.  We  get  <  gi+1(n')  <  Mj.  Thusly,  we  can 
ascertain  that  rnt  <  rnl+ 1  and  M,+ 1  <  M%.  Equivalently  restated,  we  have  a  nested  set  relation 
[mi+ 1,  Mi+ 1]  C  [rnu  M,].  As  the  iteration  i  is  allowed  to  go  to  infinity,  the  nested  set  of  possible 
values  of  gi(n)  converges  to  some  single  value.  Since  the  local  peaks  and  valleys  are  iteratively 
replaced  by  a  collasping  range  of  values  over  some  local  neighborhood,  we  aptly  name  this 
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method  the  squeeze  box  filter  (SBF).  In  our  implementation  of  the  SBF,  we  only  check  for 
Cauchy  convergence,  that  is  give  some  e  >  0  when 

OO 

X  (^(n)  “  9i+i(n))2  <  e  (17) 

n=— oo 

then  gi+\(n)  is  deemed  the  limiting  signal  or  what  we  like  to  refer  to  as  the  root  signal. 

IV.  Signal  Restoration 

A.  Additive  Noise  Restoration 

There  are  a  number  excellent  methods  to  perform  restoration  for  signals  degraded  by  the 
additive  noise  of  equation  (1)  that  are  not  constrainted  to  piecewise  constant  signals.  We  admit- 
tingly  believe  when  no  such  regularization  constraints  are  provided,  the  additive  noise  restoration 
provided  by  the  SBF  maybe  subpar  when  compared  with  methods  such  as  wavelet  thresholding, 
median  filtering,  Wiener  filtering,  OSF,  and  others.  When  constrainted  to  piecewise  constant 
sequences  the  robustness  of  the  SBF  is  exemplary  in  the  two  opposing  task  of  deceasing  the 
variation  in  intervals  where  the  ideal  signal  is  constant  and  preservation  of  edges. 

In  the  context  of  the  additive  model  given  in  equation  (1),  assuming  the  noise  has  a  stationary 
zero  mean  PDF,  to  attain  the  ideal  signal  f(n),  it  is  only  required  to  determine  the  expected 
value  of  the  detected  signal  g(n),  that  is 

E{g(n)}  =  E{p(n )  *  f(n)  +  g{n)} 

=  p{n )  *  f(n)  +  E{r](n)} 

=  Pin)*  f{n). 

Thus,  a  robust  method  to  arrive  at  E{g{n)}  yields  p(n)  *  f(n).  A  deconvolution  method  such 
as  applying  the  Weiner  filter  [13]  to  E{g{n)}  would  provide  an  estimation  for  the  ideal  signal 
f(n).  This  step  is  ignored  in  the  additive  noise  examples  in  section  V. 
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B.  Multiplicative  Noise  Restoration 

Our  proposed  restoration  method  is  not  limited  to  reconstruction  of  the  ideal  signal  polluted 
by  zero  mean  white  Gaussian  noise,  but  given  this  typical  assumption  of  the  noise,  it  will  be 
shown  that  the  blurring  caused  by  the  system  impulse  response  can  be  mathematically  accounted 
for  in  the  variance  of  the  density  function.  It  will  be  shown  that  this  translates  to  including  the 
Z2-norm  of  the  impulse  response  in  the  final  scaling  of  the  root  signal. 

The  SBF  requires  that  the  expected  value  of  the  noise  to  be  nonzero.  When  starting  with  a  zero 
mean  random  variable,  to  productively  apply  the  SBF,  the  random  variable  must  be  transformed 
so  that  the  SBF  can  be  applied  to  the  resulting  nonzero  mean  random  variable.  So  given  a  zero 
mean  normally  distributed  random  variable  X  with  variance  a\  the  PDF  is  defined  as 


Prx(x)  = 


-e 


(18) 


Ox\f^K 

where  a\  =  E{( X  —  E{X})2}.  Now,  the  random  variable  defined  by  taking  the  absolute  value 
y  =  | x |  is  a  non-Gaussian  distributed  random  variable  and  y  6  [0,  oo)  with  a  nonzero  mean, 
unless  the  variance  of  the  original  random  variable  was  zero6.  The  PDF  of  y  is  the  exponential 
distribution 

r 

PrA-(0)  if  y  =  0 


Pry  (y)  =  Pry(M)  = 

2Pr x{y)  if  0  <  y 

where  Vxx(x)  is  given  in  equation  (18).  The  expected  value  of  Y  is  determined  as 


(19) 


E{Y} 


yPxY(y)dy 


(20) 


('This  is  consider  a  noninteresting  case  and  save  for  this  footnote  is  ignore. 
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With  a  simple  substitution  (u  =  y2),  the  integral  in  equation  (20)  can  be  evaluated  and 


E{Y}  =  ax]J- 


is  easily  derived.  If  y  =  c\x\  for  some  constant  c  >  0,  then 


E{Y}  =  cax]J-. 


When  the  ideal  signal  is  equal  to  a  constant  f(n)  =  c  >  0,  the  expected  value  of  the  absolute 
value  of  the  multiplicative  model  from  equation  (2)  is 


E{\g(n)\}  =  E{\p(n)*(f(n)r](n))\} 

(  OO 

=  cE  <  p(m)r}(n  —  m) 


We  have  assumed  that  the  noise  component  is  zero  mean,  normally,  independent,  and  identically 
distribute  random  variable.  It  is  is  shown  in  [27]  that  the  weighted  sum  of  normally  distributed 


random  variables 


X  =  p{m)r]{n  —  m) 


is  a  normally  distributed  random  variable  with  a  standard  deviation  of 


<*x  —  Vn 


p2M- 


Substituting  equation  (25)  into  equation  (22)  in  the  context  of  equation  (23),  we  get 


E{\g(n)\}  =  cav  -  p2(m).  (26) 

\  7 r  z ' 

\  m=— oo 

Using  the  root  signal  of  the  SBF  to  determine  E{\g(n)\},  the  ideal  signal  f(n)  =  c  is  acquired 
by  an  appropriate  scaling  of  the  root  signal,  that  is 


/(n)  =  E{\g(n)\}  ^  P2(m) 

V  V  \m=— oo 
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V.  Examples  and  Comparative  Results 

In  this  section,  we  provide  examples  using  three  different  piecewise  constant  sequences 
corrupted  by  a  stationary  zero  mean  white  Gaussian  noise  and  a  signal  acquired  from  one  column 
of  a  Field  II  simulated  ultrasound  image  [28].  The  three  noisy  piecewise  constant  sequences  are 
either  convolved  with  the  delta  function,  a  known  finite  length  low  pass  filter  with  Z2-norm  equal 
to  one,  or  the  known  low  pass  filter  scaled  so  that  the  /2-norm  is  equal  to  three.  The  values  of 
the  unit  impulse  response  are  listed  in  table  I  and  shown  in  Fig.  5.  For  all  the  test  sequence 


p(  0) 

0.4714 

p(-l) 

0.4362 

P(l) 

pi- 2) 

0.3457 

Pi  2) 

Pi~  3) 

0.2345 

Pi 3) 

Pi~  4) 

0.1363 

Pi 4) 

Pi~  5) 

0.0678 

Pi  5) 

Pi~  6) 

0.0289 

Pi  6) 

Pi- 7) 

0.0105 

pin 

Pi- 8) 

0.0033 

Pi  8) 

Pi~  9) 

0.0009 

Pi  9) 

table  i 


Values  of  the  unit  I2-norm  impulse  response  used  in  the  test  sequences. 


77(71)  is  a  stationary  white  zero  mean  normally  distributed  random  variable,  p(n)  =  3 p(n),  and 
*  denotes  linear  convolution. 
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Fig.  5.  The  length  19  impulse  response. 


The  ideal  sequences  used  in  table  II  are 

( 


fo{n)  =  { 


1  if  0  <  n  <  99 
5  if  100  <n<  199 
1  if  200  <  n<  299, 


and 


/i  (n)  =  { 


h  (n)  =  { 


5  if  0  <  n  <  99 
1  if  100  <n<  199 
5  if  200  <n<  299, 

10  if  0  <  n  <  99 
5  if  100  <n<  199 
1  if  200  <  n<  299 
5  if  300  <  n  <  399. 


A.  A  Comparative  Qualitative  Assessment 


(28) 


(29) 


(30) 


The  original  noisy  signals  and  results  of  the  SBF  are  shown  in  Fig.s  8,  9,  10,  and  11,  resp. 
Though  exact  restoration  is  desired,  the  resulting  signal  produced  by  the  SBF  resemble  the  ideal 
piecewise  constant  signal.  It  is  evident  that  the  SBF  is  proficient  at  decreasing  the  variance  in 
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I 

gi(ri)  =  f0{n)  +  ri(n) 

II 

52(11)  =  fi(n)  +  r/(n) 

III 

gain)  =  f2(n)  +  ri(n ) 

IV 

gain)  =  \fo(n)v(n)\ 

V 

gs(n)  =  | p(n)  *  [fo(n)v(n)]\ 

VI 

gain)  =  | p{n)  *  [/o(n)»?(n)]| 

VII 

57  (n)  =  \  fi(n)r](n)\ 

VIII 

ga(n)  =  | p{n)  *  [/i(n)»?(n)]| 

IX 

g9(n)  =  | pin)  *  [/i(n)»?(n)]| 

X 

5io  (n)  =  |/i(n)?7(n)| 

XI 

5n(n)  =  1  Pin)  *  [/2(n)»?(n)]| 

XII 

5i2 in)  =  | pin)  *  [/2(n)?7(n)] 

TABLE  II 

Test  sequences. 


intervals  where  the  ideal  signal  is  a  constant  while  greatly  preserving  both  step  up  and  step  down 
edges. 

Fig.  6  is  a  Field  II  simulation  [28]  of  an  B-scan  envelop  detected  and  log  compressed 
ultrasound  image  of  five  circular  disks  of  the  same  diameter.  The  five  disks  are  representatives  of 
highly  reflective  cysts.  This  simulation  is  representative  an  image  produced  by  a  typical  clinical 
ultrasound.  Though  the  extension  of  the  SBF  to  a  two  dimensional  filter  is  straightforward,  we 
save  this  extension  for  future  publications  and  examine  only  one  column  of  the  simulated  image. 
Fig.  7(a)  we  show  the  intensity  profile  of  the  center  column  of  Fig.  6.  The  results  of  the  SBF 
and  moving  average  filter  are  in  Fig.  7(b).  The  solid  plot  of  Fig.  7(b)  is  the  root  signal  found 
by  the  SBF  and  the  dashed  plot  is  the  signal  produced  by  the  moving  average  filter.  Though 
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we  created  the  image  where  the  ideal  reflectivity  values  of  the  background  are  set  to  one  and 
the  values  within  the  disks  are  set  to  ten,  Field  II,  which  claims  a  realistic  representation  of 
an  actual  ultrasound  images,  has  rendered  a  visually  appealing  image  but  the  true  reflectivity 
values  maybe  impossible  to  recover.  Nonetheless  and  important  to  visual  acuity,  the  SBF  has 
traced  out  a  profile  that  is  large  constant  in  intervals  that  corresponds  to  the  location  occupied 
by  the  disks  while  the  large  differences  in  intensity  values  with  intervals  occcupied  by  the  disks 
and  the  background  are  preserved.  The  varying  degree  of  the  intensity  values  of  the  root  signal 
produced  by  the  SBF  in  intervals  occupied  by  the  disks  may  be  explain  by  a  unknown  spatially 
varying  PSF  in  the  original  simulated  image.  A  spatially  varying  PSF  is  a  known  phenomena 
that  occurs  with  fixed  focus  ultrasonic  imaging  systems.  Since  the  PSF  is  spatially  varying,  it 
is  reasonable  to  believe  that  the  /2-norm  of  the  impulse  response  of  the  signal  in  each  column 
of  the  simulated  image  is  varying.  Thus,  each  homogeneous  intervals  are  scaled  differently. 
This  example  serves  to  fortify  our  claim  that  for  either  additive  or  multiplicative  noisy  signals, 
which  are  representative  of  signals  produced  by  ultrasonic  imaging  systems,  the  SBF  is  robust 
in  reducing  the  variance  in  regions  where  the  signal  is  constant  while  maintaining  the  opposing 
task  of  edge  preservation. 

B.  A  Comparative  Quantitative  Assessment 

To  provide  a  quantitative  assessment  on  the  signal  restoration  performance  of  the  SBF,  the 
moving  average  filter,  the  median  filter,  Bovik  et.  al.’s  OSF  optimized  for  a  Laplacian  noise 
distribution,  Wiener  filter,  the  DWT  hard  and  soft  thresholding  methods,  we  use  a  modified  well 
known  contrast  metric  used  to  measure  contrast.  For  test  sequence  I,  II,  IV,  V,  VII,  VIII,  X,  and 
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Fig.  6.  Field  II  ultrasound  image  simulation. 


XI  to  objectively  evaluate  the  performance  of  a  restoration  algorithm  let  Q(alg)  be  defined  as 

Q(alg)  =  (31) 

01  +  ^2 

where  gai9{n )  is  the  sequence  produced  by  restoration  algorithm  ’ alg ’, 


A 4i  —  |  t  |  ^  (n)  J. 

*  nG/j 

(32) 

=  ITT  (Wn)  -  F-*)2  f 

'  ne/i 

(33) 

ii  =  [0,  99]  (J[200,  299],  and  J2  =  [100,199].  For  test  sequence  III,  VI,  and  IX,  we  use  the 
following  performance  metric 


Q(alg)  = 


(A4  i  ~~  A^)2  +  (^1  —  A/3)2  +  (^2  —  ^3)' 

0?  +  +  o’! 


(34) 


where  /it  is  defined  in  equation  (32),  of  is  defined  in  equation  (33),  I\  =  [0,99],  J2  = 
[100, 199]  1J  [300,  399],  and  J3  =  [200,299].  The  performance  metric  we  use  is  the  following 


ratio 


Q(alg) 


Qjalg) 

Q{id ) 


(35) 
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(a)  Profile  of  center  column  of 
Field  II  simulated  image. 

Fig.  7.  Intensity  profile  of  center  column  of 


(b)  SBF  scaled  (solid)  and  moving 
average  (dashed)  signals. 

;.  7,  SBF  root,  and  moving  averaged  sequences. 


where  gid(n)  =  g(n).  When  Qialg)  is  a  large  number  greater  than  one  (preferred),  then  alg 
performs  well.  If  Q(alg)  is  equal  to  one,  then  alg  does  not  provide  any  improvements  to  the 
noisy  signal.  If  Q(alg)  is  less  than  one,  then  the  tested  algorithm  alg  degrades  the  noisy  signal. 

In  our  experiment  when  performing  the  multiplicative  noise  restoration  with  Bovik’s  OSF 
optimized  Laplacian  distributed  noise,  DWT  Hard  and  Soft  thresholding,  Wiener  filtering,  we 
logarithmically  transformed  the  noisy  signal  in  an  attempt  to  cast  the  restoration  as  additive 
noise,  which  these  filters  were  design  to  overcome.  The  exponentially  transformed  values  of 
these  filters  were  then  used  in  the  quantitative  evaluation.  In  table  III  we  list  the  results  of  these 
quantitative  evaluations.  The  highest  performance  values  in  every  case  were  attained  by  the  SBF. 
This  evaluation  provide  evidence  that  the  SBF  performs  exceedingly  well. 

VI.  Conclusion 

The  additive  and  multiplicative  signal  restoration  problem  is  relevant  to  many  interesting 
applications  like  ultrasound  and  SAR  imaging.  In  this  paper,  we  constraint  the  solution  to 
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Test  sequence  — > 

I 

II 

III 

IV 

V 

VI 

VII 

VIII 

IX 

A' 

XI 

XII 

Q(MA) 

4.91 

5.3 

4.66 

9.62 

2.54 

1.98 

9.46 

2.03 

2.07 

18.0 

2.26 

1.89 

Q(Median) 

4.35 

7.79 

5.22 

4.65 

1.75 

1.41 

5.81 

1.82 

1.65 

14.0 

1.65 

1.64 

Q(OSF) 

5.24 

6.59 

5.42 

6.01 

2.0 

1.53 

6.49 

1.83 

1.64 

12.77 

1.72 

1.66 

Q(Wnr) 

1.0 

1.0 

1.0 

1.26 

0.79 

1.0 

1.18 

0.78 

0.74 

1.32 

0.98 

1.11 

Q(Hard) 

4.07 

3.19 

4.01 

4.15 

1.1 

0.0 

1.77 

1.45 

0.03 

8.64 

1.26 

0.09 

Q(Soft) 

0.60 

0.69 

0.81 

0.01 

0.02 

1.0 

0.04 

0.01 

0.82 

0.18 

0.15 

0.9 

SBF 

7.18 

7.88 

5.4 

23.01 

20.21 

20.82 

36.11 

7.58 

6.77 

71.31 

16.35 

21.56 

TABLE  III 


Q  VALUES  FOR  THE  MOVING  AVERAGE  ( MA ),  MEDIAN  (MEDIAN),  BOVIK’S  OSF  OPTIMIZED  FOR  ADDITIVE  AND 
Laplacian  distributed  noise  (OSF),  Wiener  (Wnr),  DWT  hard  thresholding  (Hard),  DWT  soft 
THRESHOLDING  (Soft),  AND  SBF  RESTORATION  METHODS. 


piecewise  constant  signals.  When  the  restoration  problem  is  viewed  in  the  stochastic  framework, 
it  can  be  analytically  shown  that  the  ideal  signal  can  be  attained  from  a  well  know  and 
easily  provable  property  of  the  expectation  operator.  The  practical  implementation  of  this  easily 
provable  property  is  problematic  and  philosophically  debateable.  We  overcome  these  problems 
by  designing  a  stochastically  driven  method  the  SBF  to  determine  the  root  signal  of  a  noisy 
signal.  In  the  additive  noise  restoration,  the  SBF  was  able  to  extract  a  exceedingly  accurate 
depiction  of  the  ideal  noise  free  signal.  The  empirical  evidence  shows  in  the  multiplicative  noise 
restoration  case  that  multiplication  with  the  analytically  derived  scaling  factor  to  the  root  signal 
yields  a  signal  that  very  well  resembles  and  very  closely  approximates  the  ideal  signal.  At  the 
very  least,  the  claim  that  the  SBF  produces  a  signal  with  reduced  variances  in  intervals  where 
the  ideal  signal  is  constant  while  the  opposing  task  of  preserving  edges  is  equally  maintained  is 
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(a)  Test  sequence  I 


(d)  Test  I:  gsBF(n)  (solid  line)  and 
gosF(n )  (dashed). 


(b)  Test  sequence  II 


(e)  Test  II:  gsBF(n)  (solid  line) 
and  gMEDiAN(n )  (dashed). 


(c)  Test  sequence  III 


(f)  Test  III:  gsBF{n)  (solid  line) 
and  gosF{n)  (dashed). 


Fig.  8.  Test  I.  Test  II.  and  Test  III  noisy,  ideal,  SBF  root,  and  next  best  Q  valued  sequences. 


subjectively  and  objectively  supported  by  the  results  of  our  experimentations. 
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0  50  100  150  200  250  300  0  50  100  150  200  250  300  0  50  100  150  200  250  300 


(d)  Test  IV:  gsBF(n)  (solid  line)  (e)  Test  V:  gsBF(n)  (solid  line)  (f)  Test  VI:  gsBF{n)  (solid  line) 

and  gMEAN(n )  (dashed).  and  gM ean (n)  (dashed).  and  g\j ean (n)  (dashed). 

Fig.  9.  Test  IV,  V,  VI  noisy,  ideal,  SBF  scaled  root,  and  moving  averaged  sequences. 
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(a)  Test  sequence  VII 


(b)  Test  sequence  VIII 


(c)  Test  sequence  IX 


(d)  Test  VII:  gsBF{n)  (solid  line) 
and  gM ean  (n)  (dashed). 


(e)  Test  VIII:  gsBF(n)  (solid  line) 
and  gM  ean  (n)  (dashed). 


(f)  Test  IX:  gsBF{n)  (solid  line) 
and  gM  ean  (n)  (dashed). 


Fig.  10.  Test  VII,  VIII,  IX  noisy,  ideal,  SBF  scaled  root,  and  moving  averaged  sequences. 
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(a)  Test  sequence  X 


(b)  Test  sequence  XI 


(c)  Test  sequence  XII 


(d)  Test  X:  gsBF{n)  (solid  line)  (e)  Test  XI:  gsBF{n)  (solid  line)  (f)  Test  XII:  gsBF(n)  (solid  line) 

and  gMEAN(n )  (dashed).  and  gMEAN(n )  (dashed).  and  gM ean  (ft)  (dashed). 

Fig.  11.  Test  X.  XI,  XII  noisy,  ideal,  SBF  scaled  root,  moving  averaged  sequences. 
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Abstract — Reverberation  and  multi-path  reflection  artifacts 
are  a  common  problem  in  ultrasound  imaging.  We  propose  a 
novel  method  to  remove  these  artifacts.  Regions  adversely  af¬ 
fected  by  these  artifacts  are  replaced  with  textures  that  resemble 
the  underlying  object(s),  which  were  originally  obscured.  Our 
proposed  method  incorporates  optimally  soft  thresholding  the 
2D  discrete  wavelet  transform  of  the  artifact  regions  to  produce 
a  near  optimal  estimate  of  the  reflectivity  values  due  only  to 
the  reverberation  and  multi-path  reflection  artifacts.  Simply 
subtracting  this  estimate  from  the  original  reflectivity  values, 
we  attain  a  near  optimal  estimate  of  the  artifact  free  reflectivity 
values.  We  provide  B  mode  images  to  substantiate  the  benefits 
of  this  method  in  producing  a  more  useful  and  visually  pleasing 
image. 

I.  Introduction 

The  use  of  ultrasound  to  provide  a  useful  imaging  tool  is 
preferred  because  of  the  noninvasive  and  nonionizing  nature 
of  this  realtime  imaging  system  that  can  provide  functional 
information  like  motion  and  velocities.  Although  the  images 
and  videos  produced  by  a  ultrasound  system  are  generally 
beneficial,  its  usefulness  as  a  tool  to  aid  medical  diagnosis  can 
be  improved  by  removing  the  various  artifacts  that  may  arise. 
Two  types  of  ultrasound  artifacts  are  caused  by  reverberation 
and  muli-path  reflection  of  the  ultrasonic  sound  wave  as  it 
travels  across  or  around  highly  reflective  objects  or  interfaces. 
The  description  of  these  artifacts  given  in  [1]  are  paraphased 
in  the  following: 

•  Reverberation  artifacts  occur  when  reflected  ultrasound 
energy  is  reflected  back  and  forth  between  two  closely 
spaced  interfaces  during  signal  acquistion  and  prior  to 
the  next  transmitted  pulse. 

•  Artifacts  due  to  multi-path  reflection  occur  when  the 
ultrasound  beam  is  nonperpendicularly  reflected  or  re¬ 
fracted  off  a  highly  reflective  surface  and  subsequently 
detected  at  the  transducer. 

Since  the  range1  is  directly  related  to  time  in  ultrasound 
imaging,  reverberation  artifacts  can  be  seen  as  multiple  equal 
spaced  objects  with  amplitudes  decreasing  as  depth  increases. 
Also,  the  paths  of  the  nonperpendicularly  reflected  beams 
are  longer  than  the  paths  of  perpendicularly  reflected  beams. 
Consequently,  highly  reflective  objects  may  reappear  further 
away  from  the  transducer  and  the  multi-path  refection  artifacts 

'The  distance  from  the  transducer  or  depth. 


are  seen  as  misplaced  objects.  In  ultrasound  images  of  organs, 
muscles,  or  tissues  that  are  in  close  proximity  of  highly 
reflective  structures  such  as  bones  or  tissue/gas  interfaces, 
the  undesirable  occurrences  of  reverberation  and  multi-path 
reflection  artifacts  are  common. 

An  example  application  to  show  the  importance  of  the 
removal  of  these  artifacts  is  improving  the  ultrasound  imaging 
system  so  that  a  3D  rendering  of  muscles  and  tissues  within 
the  heart  can  be  better  observed  and  important  information 
such  as  the  volume  and  surface  area  of  the  heart  can  be 
accurately  assessed.  In  principle,  by  specifically  removing  the 
reverberation  and  multi-path  reflection  artifacts  caused  by  the 
highly  reflective  nature  of  the  ribs  we  offer  improvements  to 
the  B  mode  or  M  mode  ultrasound  images  so  that  the  motions 
and  borders  of  the  myocardium2  of  the  heart  can  accurately 
be  evaluated  by  an  automated  algorithm,  a  physician,  or  both. 
Additionally,  a  3D  rendering  that  depicts  the  surface  of  the 
heart  and  volume  information  using  these  2D  images  is  not 
adversely  skewed  by  these  common  artifacts. 

An  appropriately  placed  short  axis  scans  of  the  left  ventricle 
will  show  the  myocardium  contracting  and  expanding  as  the 
heart  repeats  its  end  diastole  to  end  systole  cycle.  In  acquiring 
a  short  axis  view  of  the  heart,  the  ultrasonic  sound  wave  must 
traverse  through  or  around  the  ribs.  The  highly  reflective  nature 
of  the  ribs  can  cause  reverberation  and/or  multi-path  artifacts 
to  appear  within  the  region  of  interest  that  is  over  or  near 
the  mycardium.  The  presence  of  these  artifacts  can  obscure 
parts  of  the  myocardium  and  makes  determining  its  motions 
and  borders  problematic.  We  attempt  to  illustrate  the  problem 
in  Fig.  2(a)  where  a  typical  short  axis  B  mode  ultrasound 
image  of  the  left  ventricle  of  a  mouse  heart  is  shown.  In 
Fig.  2(a)  the  approximated  location  of  the  inner  left  ventricle 
myocardial  border  is  shown  by  the  ellipse.  The  reverberation 
or  multi-path  reflection  artifacts  most  probably  caused  by  the 
ribs  are  are  seen  as  bright  objects.  We  have  highlighted  these 
artifacts  with  a  bold  “A”  in  Fig.  2(a).  A  B  mode  image 
that  isolates  these  artifacts  is  shown  in  Fig.  2(b).  It  can  be 
observed  that  the  leftmost  section  of  the  myocardium  are  made 
ambiguous  by  these  artifacts.  To  accurately  assess  the  health 
of  the  myocardium  or  other  obscured  muscle(s)/tissue(s),  it  is 

-The  middle  and  thickest  layer  of  the  heart  wall,  composed  of  cardiac 
muscles  [2], 


necessary  not  only  to  remove  these  artifacts,  but  the  reflectivity 
values  of  the  muscle(s)/tissue(s)  obscured  by  these  artifacts 
should  be  accurately  preserved. 

II.  Background 

A.  The  Discrete  Wavelet  Transform 

In  [3],  Mallat’s  multiresolution  wavelet  transform  provided 
a  two  dimensional  (2D)  subband  decomposition  of  an  image 
that  allows  perfect  reconstruction.  Although  our  proposed 
artifact  removal  method  strives  to  remove  certain  components 
of  the  original  image,  perfect  reconstruction  is  still  preferred 
since  no  visual  information  loss  due  to  the  transform  is 
guaranteed.  Since  one  of  our  proposed  goal  is  to  render  a 
more  visually  pleasing  reverberation  and  multi-path  reflection 
artifact  free  image,  our  DWT  construction  emulates  the  forty- 
three  channel  Gabor  filter  bank  of  [4].  The  2D  filter  bank 
composed  of  2D  Gabor  filters  in  [4]  was  motivated  by  physio- 
psychological  experimental  evidence  that  the  early  stages  of  a 
biological  vision  system  are  well  represented  by  conjointly 
well  spatial  and  frequency  localized  bank  of  Gabor  filters 
where  the  number  of  channels  is  in  the  forties.  Another  key 
characterization  of  this  Gabor  filter  bank  construction  is  that 
the  magnitude  response  of  the  filters  becomes  wider  and  their 
magnitudes  decrease  as  the  center  hortizontal  and/or  vertical 
frequencies  increase.  With  these  characterizations  in  mind  we 
emulate  the  Gabor  filter  bank  in  [4]  with  a  2D  forty  channel 
well  conjointly  localized  perfect  reconstruction  DWT  filter 
bank.  The  forty  channel  2D  DWT  decomposition  is  shown  in 
Fig.  ??  and  for  comparison  purpose  the  typical  thirteen  four 
level  DWT  decomposition  is  shown  in  Fig.  ??.  Evidence  is 
provide  in  [5]  that  shows  the  multi-level  DWT  filter  bank  is 
conjointly  well  localized  when  the  Coifman  quadrature  mirror 
filter  bank  (QMF),  which  maximizes  the  number  vanishing 
moment  for  a  given  support  width,  is  used.  A  description  of 
the  Coifman  QMF  used  in  our  2D  forty  channel  DWT  filter 
bank  can  be  found  in  [6]. 

B.  Hard  and  Soft  Thresholding  Methods 

For  signals  corrupted  by  additive  white  Gaussian  noise  such 
as  in  equation  (1),  a  typical  denoising  method  to  recover  the 
noise  free  signal  x(n)  is  to  apply  a  hard  or  soft  thresholding 
operator  to  the  DWT  coefficients.  These  methods  are  com¬ 
monly  referred  to  as  wavelet  shrinkage.  The  hard  and  soft 
thresholding  operation  are  given  in  equations  (2)  and  (3),  resp. 
where  w  =  DWT{ y}. 

y(n)  =  x(n)  +  rj(n)  (1) 

The  noise  free  signal  is  reconstructed  from  the  threshold 
wavelet  coefficents  y  =  IDWT{ w}  «  x  where  IDWT{-} 
means  inverse  DWT  (IDWT). 

f. _  /  0  if  lU;(n)l  <  ^ hard,  /o) 

\  w(n)  otherwise. 

!w(n)  -  Xsoft  if  w(n)  >  X soft 

w(n )  +  A soft  if  w(n)  <  -Xsoft  (3) 

0  otherwise. 


The  robustness  of  the  wavelet  shrinkage  method  is  depen¬ 
dant  on  the  choice  of  threshold(s).  Some  various  methods 
and  criteria  to  consider  in  choosing  a  global  or  level  depen¬ 
dant  hard  or  soft  threshold(s)  based  on  a  priori  known  or 
estimated  noise  statistics,  signal  length,  mean  squared  error, 
and  smoothness  are  described  in  [7].  In  [7],  [8],  and  [9],  the 
wavelet  shrinkage  denoising  of  the  signal  y(n)  using  hard  or 
soft  thresholding  of  the  DWT  coefficient  with  a  global  soft 

threshold  value  _ 

Xsoft  =  a  \/  2  ln(iV) ,  (4) 

or  some  variant  of  this  threshold  where  cr  is  the  standard 
deviation  of  the  noise  and  N  is  length  of  the  signal,  is  shown 
to  be  optimal  or  near  optimal  in  the  sense  that  the  mean 
squared  error  is  minimized  (MMSE)  or  the  maximum  of  the 
mean  squared  error  is  minimized  (minimax).  Additionally  the 
reconstructed  signal  of  the  soft  threshold  wavelet  shrinkage 
method,  Donoho  in  [8]  claim  is  nearly  as  smooth  as  the 
original  signal  where  the  smoothness  of  the  reconstructed 
signal  was  determined  from  a  “wide  range”  of  smoothness 
metric  [8]. 

As  a  stringent  motivation  in  our  ultrasound  artifact  removal 
algorithm,  we  desire  mean  squared  error  optimality  or  at  least 
near  optimality  in  the  sense  of  MMSE  or  minimax.  To  prevent 
processing  artifacts  such  as  unwarranted  oscillations  (ringing) 
or  discontinuities,  we  adhere  to  the  constraint  that  the  recon¬ 
structed  signal/image  should  be  equivalently  as  smoothness  as 
the  original  signal/image.  We  will  take  advantage  of  the  work 
accomplished  in  [7],  [8],  and  [9]  by  using  the  soft  threshold 
wavelet  shrinkage  method  with  the  threshold  value  defined  as 
in  equation  (4)  in  our  proposed  algorithm. 

III.  The  Artifact  Removal  Method 

In  [10]  Bjaerum  and  Torp  proposed  an  additive  model  of  the 
complex  demodulated  Doppler  signal  to  remove  clutter,  which 
are  objects  that  do  not  move  or  move  slowly.  In  modeling  the 
complex  demodulated  IQ  data,  we  adopted  a  similar  additive 
model.  The  model  of  the  IQ  data  that  we  have  adopted  is  the 
sum  of  three  complex  value  components 

IQ(n,  to)  =  A(n,  to)  +  T(n,  in)  +  77(71,  to)  (5) 

where  A(n,  to)  is  the  complex  value  due  to  the  artifacts, 
T(n,  to)  is  the  complex  reflectivity  value  due  to  the  underlying 
muscle(s)/tissue(s)  obscured  by  A(n,  to),  and  77(71,  m)  is  com¬ 
plex  valued  white  noise  where  the  real  and  the  imaginary  parts 
of  77(71,  to)  are  zero  mean  Gaussian  distributed.  It  is  reasonable 
to  assume  that  the  complex  values  due  to  the  artifacts  and  the 
reflectivity  values  due  to  the  underlying  muscle(s)/tissue(s)  are 
independent  of  each  other. 

The  B  mode  data  is  attained  by  log  compressing  the  L2- 
norm  of  the  interpolated  IQ  data.  The  data  flow  of  the  IQ 
data  to  B  mode  data  is  shown  in  Fig.  1.  Although  the  B 
mode  data  is  used  to  produce  a  visually  meaningful  image  for 
medical  diagnosis,  processing  the  IQ  data  is  better  suited  for 
our  current  application  of  removing  reverberation  and  multi- 
path  reflection  artifacts.  It  is  more  advantageous  to  process  the 
IQ  data  instead  of  the  B  mode  data  for  the  following  reasons: 


B  mode 


and 


IQ- 


HI 

Interp(-) 

log(-) 

Fig.  1.  IQ  data  to  B  mode  image. 


•  The  values  of  the  IQ  data  encompass  a  greater  dynamic 
range  than  the  B  mode  data.  Thus,  samples  that  are 
greatly  affected  by  artifacts  are  more  easily  distinguish¬ 
able. 

•  The  number  of  samples  of  IQ  data  is  substantially  less 
than  the  number  of  B  mode  data.  Always  important  to 
real  time  or  near  real  time  algorithms  is  to  process  the 
fewest  number  of  samples. 

•  The  signal  to  noise  ratio  (SNR)  of  the  IQ  data  is  greater 
than  the  SNR  of  the  B  mode  data  where  SNR  is  meant  as 
the  ratio  of  the  power  due  to  the  underlying  muscle/tissue 
reflectivity  values  T(n,m)  and  the  power  due  to  white 
noise  77(71,  m).  The  reduction  in  SNR  of  the  B  mode  data 
is  due  to  log  compression,  which  is  necessary  to  render  a 
image  that  is  within  the  dynamic  range  of  human  vision. 

•  The  real  and  imaginary  part  of  T(n,m )  are  Gaussian 
distributed.  This  detail  is  required  for  the  optimality  or 
near  optimally  of  the  soft  thresholding  wavelet  shrinkage 
method  using  the  threshold  given  in  equation  4.  We  will 
elaborate  and  establish  this  detail  later  in  this  section. 


From  the  additive  model  we  propose  in  equation  (5),  if 
a  sample  is  adversely  affected  by  the  artifact  component, 
then  the  L2-norm  of  the  artifact  component  is  substantially 
greater  than  the  L2-norm  of  the  muscle/tissue  component,  that 
is  ||A(n, m)\\  ||T(n, m)\\.  An  estimate  of  samples  where 

artifacts  are  prominent  can  be  determined  by  thresholding  the 
L2-norm  of  the  IQ  data 


IQ(n,  m) 


IQ(n,m )  if  \\IQ(n,m)\\  >  A 
0  otherwise. 


(6) 


The  estimated  artifact  dominanted  complex  data  IQ  contains 
the  complex  reflectivity  values  of  the  underlying  muscle  and/or 
tissue,  which  are  obscured  by  the  artifact(s).  To  provide  an 
estimate  of  the  reflectivity  values  of  the  underlying  muscle 
and/or  tissue,  we  consider  the  complex  reflectivity  values  of 
the  underlying  muscle  or  tissue  T(n,  m)  of  the  estimated  arti¬ 
fact  dominanted  complex  date  IQ(n,  m )  as  unwanted  additive 
noise,  that  is 

IQ(n ,  m)  =  A(n,  in)  +  T(n,  m)  +  77(71,  m) .  (7) 

S,  y  ^ 

noise 

Applying  a  wavelet  shrinkage  algorithm  to  the  real  and  imag¬ 
inary  parts  of  IQ  to  remove  T  and  rj ,  we  attain  an  estimate 
of  the  complex  values  due  only  to  the  artifact  A.  More 
precisely,  let  W reai  =  DWT{real{ IQ}}  and  W irnag  = 
DWT{imag{IQ}}  be  the  2D  DWT  of  the  real  and  imaginary 
image  of  IQ,  resp.  The  soft  thresholding  of  Wrea/  and  W imag 
with  thresholds 


A 


so  ft, real  tT real 


^2  IndTQI) 


^so  ft,  imag  ®imag  y  21n(|IQ|),  (9) 

resp.,  are  denoted  as  W reai  and  W imag ,  resp.  The  term  <jreai 
is  the  standard  deviation  of  elements  in  IQreai  where 

IQreai  =  {real{IQ(n,m)}  :  \\IQ(n,m)\\  <  A}  (10) 

with  the  same  A  used  in  equation  (6).  Likewise  aimag  is  the 
standard  deviation  of  IQimag  where 

IQimag  =  {imag{IQ(n,m)}  :  \\IQ(n,m)\\  <  A}  (11) 

with  the  same  A  used  in  equation  (6).  The  term  |IQ|  in  equa¬ 
tions  (8)  and  (9)  is  simply  the  number  of  samples  where  the 
L2-norm  of  the  IQ  data  is  greater  than  the  A  of  equation  (6). 

An  estimate  of  the  artifacts  in  equation  (7)  is  determine  as 

A  =  IDWT{  Wreal}  +  jIDWT{  Wimag}.  (12) 

A  artifact  free  image  IQ  is  attained  by  subtracting  the  esti¬ 
mated  artifact  values  of  equation  (12)  from  the  original  IQ 
data,  that  is 

IQ  =  IQ  A.  (13) 

If  the  component  T(n,m )  in  the  additive  model  given  in 
equation  (7)  can  be  shown  to  be  Gaussian  distributed,  then 
using  the  soft  threshold  of  equation  (4),  which  we  estimate  in 
equations  (8)  and  (9),  would  yield  a  near  optimal  estimate  of 
the  artifact  values  A(n,m).  If  the  component  T(n,m )  were 
statistically  independent  with  respect  to  the  sample  indices 
(n,  m)  £  Z  x  Z,  then  optimally  in  the  sense  describe  by 
Donoho  et  al.  in  [7],  [8],  and  [9]  could  justifiably  be  claimed3. 

To  establish  that  T(n,m)  is  Gaussian  distributed,  we  con¬ 
sider  the  basic  characterization  of  these  reflectivity  values. 
In  [11]  and  [12],  Goodman  characterize  the  reflectivity  values 
produced  by  ultrasound,  coherent  optical  laser,  and  synthetic 
aperture  radar  as  a  sum  of  complex  random  phasors 

OO 

T(n,  m)  =  ^  a,;(n,  m)eJv^n,rra^  (14) 

i=0 

where  cii(n,m )  and  <£>j(n,m)  are  independent  with  respect 
to  the  variable  t  £  Z  and  with  each  other4.  This  specular 
phenomena  occurs  naturally  and  normally  with  these  and  other 
imaging  systems.  The  speckle  characterization  is  due  to  the 
“roughness”  of  the  object  being  imaged  with  respect  to  the 
wavelength  of  the  transmitted  sound,  light,  or  electro-magnetic 
wave.  Since  T(n,m)  is  the  infinite  sum  of  independent 
variables,  the  Central  Limit  Theorem  implies  that  T(n,m .)  is 
Gaussian  distributed. 

3  Since  the  reflectivity  values  due  to  the  underlying  muscle/tissue  are  not 
statistically  independent,  the  arguments  for  MMSE  or  minimax  optimality 
made  by  Dohono  et  al.  cannot  be  invoked. 

4  For  fixed  n,m,i  G  the  conditional  probability  of  ai(n,m)  given 
( pi(n,m )  is  equal  to  the  of  the  probablity  of  aj(n,m)  and  likewise  for  the 
conditional  probability  of  (pi(n ,  m)  given  ai(n,  m)  is  equal  to  the  probability 
of  (pi(n ,  m). 


(8) 


IV.  Results 

A.  Image  Artifact  Removal 

In  Fig.  2  we  show  the  B  mode  images  of  our  proposed 
reverberation  and  multi-path  artifact  removal  algorithm.  The 
images  shown  in  Fig.  2  are  rendered  in  the  same  dynamic 
range  so  that  an  accurate  representation  of  our  results  can 
be  displayed.  Although  our  processing  is  performed  on  the 
IQ  data,  our  end  result  is  to  improve  the  quality  of  the  B 
mode  image  and  we  use  this  mode  to  display  our  results. 
The  B  mode  image  in  Fig.  2(a)  is  a  short  axis  view  of  a 
mouse  heart  where  the  bold  capital  “A”  signifies  the  presence 
of  an  artifact  and  the  ellipse  represents  the  approximate 
location  of  the  muscle/tissue  of  interest.  In  Fig.  2(b)  is  the 
B  mode  image  of  IQ(n,m .)  as  defined  in  equation  (6)  using 
A  =  400.  The  image  in  Fig.  2(b)  contains  both  prominent 
artifacts  and  muscle(s)/tissue(s),  which  are  obscured.  The  B 
mode  image  of  the  near  optimally  estimated  artifact  only 
reflectivity  values  A  of  equation  (12)  is  shown  in  Fig.  2(c). 
When  compared  with  the  image  in  Fig.  2(b),  the  objects  and 
features  in  the  B  mode  image  of  Fig.  2(c)  are  subtly  and 
smoothly  diminished.  Importantly,  it  is  evident  from  Fig.  2(c) 
that  no  new  discontinuity,  ripples,  blips,  or  oscilliations  are 
introduced  by  our  proposed  processing.  The  B  mode  image 
of  the  approximated  artifact  free  image  produced  by  our 
algorithm  IQ(n,m )  is  shown  in  Fig.  2(d).  We  have  yet  to 
make  any  objective  and  quantifiable  measurements  on  the 
robustness  of  our  proposed  artifact  removal  algorithm.  From  a 
subjective  evaluation  of  Fig.  2(d),  we  highlight  the  following 
improvements  and  observations: 

1)  The  reverberation  and  multi-path  artifacts,  which  are 
evident  in  Fig.  2(a),  are  removed. 

2)  The  artifact  regions  have  been  replace  by  textures  that 
are  homogeneous  with  textures  from  neighboring  re¬ 
gions  that  were  not  adversely  affected  by  these  artifacts. 

3)  Artifact  free  regions  are  not  diminished  by  our  proposed 
algorithm. 

4)  No  processing  artifacts  are  visibly  evident. 

5)  Objects  of  interest  are  not  morphed,  warped,  skewed,  or 
disfigured  in  anyway. 

6)  The  image  in  Fig.  2(d)  is  not  displayed  in  the  full 
contrast  range  of  the  human  visual  system.  A  contrast 
enhancement  would  render  a  more  visually  pleasing 
image. 

B.  3D  Surface  Rendering  and  Volume  Estimation 

We  implement  a  3D  extenstion  of  the  2D  gradient  vector 
flow  (GVF)  active  contour  (also  known  as  a  snake)  (cite  Xu 
and  Prince)  to  preform  a  3D  rendering  from  2D  slices.  The 
volume  of  the  3D  rendering  is  calculated  as  the  volume  of 
convex  hull  of  the  final  3D  GVF  snake.  In  Fig.  we  shows  the 
3D  contour  using  unprocessed  slices  and  the  volume  of  this 
contour  is  estimated  at  (some  number).  In  Fig.  we  show  the  3D 
rendering  provided  by  the  3D  GVF  snake  and  using  artifact 
removed  slices.  The  volume  enclosed  by  the  3D  contour  in 
Fig.  is  (some  more  accurate  number).  This  result  is  more 


accurate  to  the  true  volume.  Although  our  evaluation  is  far 
from  complete,  this  one  example  shows  the  artifact  removal  is 
necessary  for  accurate  3D  imaging  and  our  proposed  artifact 
removal  method  is  promising. 

V.  Conclusion 

Ultrasound  artifacts  due  to  reverberations  or  multi-path 
reflections  are  expected  when  acquiring  images  of  organs, 
muscles,  tissues,  etc.  require  the  ultrasonic  sound  wave  to  tra¬ 
verse  through  or  around  highly  reflective  objects  such  as  bones 
or  various  interfaces.  We  present  a  wavelet  transform  method 
to  replace  these  artifacts  with  a  near  optimal  estimate  of  the 
underlying  objects,  which  are  obscured  by  these  artifacts.  For 
several  compelling  reasons,  our  processing  is  performed  fully 
to  the  complex  IQ  data.  Using  the  resulting  B  mode  data 
produced  by  the  processed  IQ  data,  we  show  that  the  regions 
adversely  affected  by  artifacts  are  replaced  with  textures  that 
are  homogeneous  with  textures  from  surrounding  regions  not 
adversely  affected  by  these  artifacts.  Additionally,  it  can  be 
observed  that  our  proposed  artifact  removal  algorithm  is  not 
detrimental  to  artifact  free  regions  and  no  processing  artifacts 
are  introduced. 
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(a)  The  original  image  with  artifacts  (A)  and  approximate 
location  of  myocardium  (ellipse). 


(b)  Unprocessed  artifact  image  IQ  using  A  =  400 


(c)  Estimated  artifact  only  image  A.  (d)  Artifact  free  image  IQ(n,m) 


Fig.  2.  B  mode  version  of  original,  unprocessed  artifact,  processed  artifact,  and  artifact  free  images. 
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ABSTRACT 

A  novel  stochastically  driven  filtering  method  to  despeckle  B 
mode  ultrasound  images  is  presented.  This  method  is  motivated  by 
viewing  the  pixel  values  as  a  stochastic  process  and  removing  out¬ 
liers,  where  outliers  are  defined  by  local  extrema.  These  outliers 
are  removed  by  local  averaging.  This  produces  another  image  with 
new  outliers  (local  extrema)  and  the  process  is  iteratively  repeated. 
With  each  iteration  homogeneous  regions  become  smoother  while 
edges  that  defined  these  regions  remain  preserved.  To  evaluate  the 
performance  of  our  proposed  method  in  satisfying  these  two  op¬ 
posing  goals  we  develop  a  modified  Fisher  discriminant  contrast 
metric.  Larger  values  of  this  metric  indicate  better  performance  in 
reducing  each  intraregion  or  intraclass  variance  and  increasing  the 
difference  of  interregion  or  interclass  means. 

1.  INTRODUCTION 

In  applications  where  speckle  precludes  successful  image  analysis 
and  removal  is  desired,  speckle  is  considered  noise  and  its  removal 
as  an  image  restoration  problem.  Speckle  is  a  common  phenom¬ 
ena  found  in  many  imaging  modalities  such  as  optical  laser,  syn¬ 
thetic  aperture  radar,  and  ultrasound.  Many  despeckling  methods 
have  been  proposed  with  improving  these  imaging  modalities  in 
mind.  The  despeckling  success  of  the  various  proposed  algorithms 
are  usually  subjectively  assessed.  A  fair  quantitative  evaluation  in 
many  cases  is  avoided. 

Our  concern  is  despeckling  B  mode  ultrasound  images  to  aid 
medical  diagnosis.  Wagner  et  at.  in  [1,  2]  gives  a  description  of 
the  statistical  characteristics  of  B  mode  ultrasound  speckle.  The 
assumptions  of  the  cause  and  basic  characterizations  of  ultrasound 
speckle  coincide  with  Goodman's  assumptions  and  resulting  sta¬ 
tistical  characterizations  of  speckle  caused  by  a  coherent  laser 
in  [3,  4],  Although  these  characterizations  of  speckle  are  insight¬ 
ful,  an  image  model  like  the  additive  and  multiplicative  ones  in 
equations  (1)  and  (2),  resp.,  is  not  offerred.  Two  typical  image 
models  are 

J(n,m)  =  I{n,m)  +  rj(n,m)  (1) 

for  the  additive  noise  case  and 

J(n,m)  =  I(n,m)rj(n,m)  (2) 

for  the  multiplicative  noise  case  where  r}(n ,  m)  is  the  noise  com¬ 
ponent.  Most  image  restoration  methods  are  specific  to  the  prob¬ 
lem  being  addressed.  The  restoration  problem  or  motivation  for  a 
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reasonable  solution  are  generally  constrained  by  the  modeling  of 
the  noisy  image. 

We  give  a  brief  overview  of  various  despeckling  algorithms 
proposed  for  different  imaging  modalities  and  applications.  A 
novel  stochastically  driven  method  to  remove  or  reduce  ultrasound 
speckle  is  presented.  Lastly,  we  provide  a  quantitative  assessment 
of  the  results  of  our  proposed  method  and  those  of  other  despeck¬ 
ling  methods. 

2.  VARIOUS  DESPECKLING  METHODS 

Nagao  and  Matsuyama  in  [5]  proposed  a  recursive  edge  preserving 
smoothing  algorithm.  One  iteration  of  this  algorithm  replaces  each 
pixel  value  with  the  mean  of  some  segment  wu  originating  from 
■J(n,  rn )  and  the  variance  of  wk  is  the  minimum  variance  attained 
from  a  set  of  variances  in  all  eight  directions.  Precisely, 

J (n,  m)  =  E(wk),  (3) 

where  £(•)  is  the  expected  value  operator, 

var(wk)  =  min {var(w,)  |  fori  =  0, 1,  2, . . . ,  7} 

and  wo,  Wi,  W2,  ■  ■  ■ ,  w?  are  equal  length  segments,  which  origi¬ 
nate  from  J(n,  m)  and  span  all  eight  directions.  This  recursion  is 
allowed  to  continue  until  convergence  in  most  of  the  pixel  values 
is  established. 

Lee  in  [6]  proposed  methods  to  contrast  enhance  an  image 
and  to  restore  an  image  corrupted  by  noise  as  define  in  equa¬ 
tions  (1)  or  (2).  The  method  in  [6]  proposes  to  adaptively  locally 
smooth  in  homogeneous  regions  while  regions  containing  edges  or 
salient  textures  are  preserved.  This  algorithm  adaptively  chooses 
a  weighting  factor  k  between  zero  and  one  so  that  the  new  pixel 
value  J(n,  m)  is  set  at 

■J(n,  m)  =  p  +  k(J(n,  m)  —  p)  (4) 

where  p  is  the  mean  in  some  window.  The  weighting  factor  k  is 
adaptively  determined  as 


where  a2  is  the  local  variance  in  the  some  window.  The  noise 
variance  is  denoted  as  a 2  and  must  be  known  a  priori.  When 
ff2  at,  A  0,  the  gain  parameter  is  approximately  one,  in  which 
case  the  filter  in  equation  (4)  is  the  identity  filter.  If  the  local  vari¬ 
ance  a2  is  greater  than  but  nearly  equal  to  the  global  noise  variance 


a2,  then  the  filter  specified  by  equation  (4)  is  a  local  averaging  low 
pass  filter. 

When  the  image  is  degraded  by  multiplicative  noise  as  in  equa¬ 
tion  (2),  then  Lee  in  [6]  recasts  the  image  as  an  additive  noise 
model  of  the  form 

J{n,  m)  =  AI(n,  m)  +  Brj(n,  m)  +  C  (6) 


where  w(n,m')  =  1  and  W  is  some  odd  dimensional 

n'  ,m' 

square  window. 

The  adaptive  weighted  median  filter  (AWMF)  of  Loupas  et  al. 
in  [9]  models  the  ultrasound  image  as 

■J{n,  m)  =  I{n,  m)  +  \/l(n ,  m)q(n,  m).  (16) 


where  A,  B,C  £  R  are  chosen  so  that  the  mean  squared  error  of 
the  multiplicative  model  from  equation  (2)  and  the  additive  model 
in  equation  (6)  is  minimized.  Recasting  as  an  image  with  additive 
noise,  the  adaptative  filter  defined  in  equation  (4)  can  be  used  and 
the  gain  parameter  k  is  redefined  as 


pTj  Q 

njtf  +  h2Q 


(7) 


where  //  /  is  the  approximated  local  mean  of  the  ideal  image 
I  (n,  m )  determined  as 

Hi  =  —  (8) 

Hv 

within  some  fixed  window  and  p  is  the  local  mean  of  J(n,  m). 
The  variable  Q  is  determined  as 


Q 


a2  +  H2  2 
<*;;  +  Hv  111 


(9) 


where  a2  is  the  local  variance  of  J(n,  m)  within  some  window. 

The  Frost  filter  given  in  [7]  consists  of  determining  a  filter 
f ( n  .  in'),  so  that  in  a  local  homogeneous  region  of  J{n,  rn )  and 
in  the  presence  of  white  noise,  the  expection  criterion 


E[(I(n,  m)  —  I(n,  m))1]  (10) 


is  minimized.  The  term  I(n,m  )  is  a  windowed  weighted  sum 
about  some  constrained  ideal  image  I(n,m.) 


I(n,  m)  =  E  f(n  ,  m')I(n  +  n  ,  m  +  m)  (11) 

n'  ,m'  G  W 


where  W  is  some  odd  dimensional  window.  They  derived  the  filter 
that  minimizes  equation  (10)  as 

f(n,m')  =  Kae~°  (12) 

where 


K  is  some  normalizing  constant,  a  is  an  region  dependant  con¬ 
stant,  fiv  is  the  mean  of  the  noise,  av  is  the  standard  deviation  of 
the  noise,  fi  and  a  are  the  local  mean  and  local  standard  devia¬ 
tion  of  J(n,  m),  resp.  The  function  r  :  7L  x  7L  — ►  R  is  some 
symmetric  function  like  r(n,  rn)  =  y/n2  +  m,2. 

The  linear  minimum  mean  squared  error  estimate  of  the  ideal 
noise  free  image  given  by  Kuan  et  al.  in  [8]  is  defined  as 

2  2 

.  LO  —  (Jr,  .  ,  ,  , 

/(n,  m)  =  /j,  -\ - 2-(  J(n ,  m)  -  fi)  (14) 

LOz 

where  //  is  the  local  variance  in  some  window  about  J(n,  m)  and 
Ld2  is  the  local  variance  in  some  weighted  window  about  J(n,  m). 
The  local  weighted  variance  lo2  is  defined  as 

u2  =  N}M,  w(n  ,m')(J(n  +  n  ,m.  +  m')-p.)2  (15) 

n'  ,m'  G  W 


From  the  image  model  given  in  equation  (16)  and  when  the  ideal 
image  is  equal  to  a  constant  in  some  local  neighborhood,  then 

a  =  cav  (17) 

where  a2  and  a %  are  the  local  variance  of  the  observed  image  and 
the  variance  of  the  noise,  resp.  The  AWMF  is  defined  as 

J(n,  m)  =  median (  J{n  ,  m'), . .  . ,  J{n  ,  m')  )  (18) 

w  (n1  ,mr ) 

for  n  ,  m  in  some  window  about  (n,  m)  and 

/  '  ,  f  m  m  co2 \/n'-  +  m'2  } 

w(n  ,m)  =  round+  <  ic(U,  U) - >  (19) 

l  /'  J 


where  round+{-}  means  round  to  the  nearest  non-negative  inte¬ 
ger,  c  is  some  constant,  p  is  the  the  local  mean,  and  a2  is  the  local 
variance.  The  constant  c  and  the  window  term  u;(0,  0)  determines 
the  AWMF's  ability  to  preserve  edges. 

Yu  and  Acton  in  [10]  proposed  the  speckle  reducing 
anisotropic  diffusion  (SRAD)  method.  The  SRAD  algorithm  it¬ 
eratively  filters  a  nonzero  valued  image  I{x,  y,  0)  =  I(x,  y)  ac¬ 
cording  to 

=  div[c{q)X7I{x,y,t)]  (20) 


and 


dl(x,  y,  t)  =  Q 

where  n  is  the  outward  normal  vector  to  the  border  of  fi. 
diffusion  coefficient  c(-)  is  defined  either  as  the  quotient 


(21) 

The 


c(q(x,y;t)) 


qUt)  +  q6(t) 

q£(t)  +  q2{x,y,t) 


(22) 


or  as  the  exponential  function 


c(q(x,y,t))  =  exp 


( qo(t)  ~  q2(x,y,t)\ 

V  q40(t)  +  q2(t)  )■ 


(23) 


In  both  equations  (22)  and  (23),  if  q(x,y,t)  «  qo(t),  then 
c(q(x,y,t))  ~  1  and  equation  (20)  is  a  local  smoothing  opera¬ 
tor.  If  q(x,y,t)  qo(t),  then  the  diffusion  coefficient  is  very 
small  and  smoothing  in  a  local  region  around  (x,  y)  is  suppressed. 
When  at  time  t,  if  (x,  y)  resides  in  a  homogeneous  region,  then 
smoothing  can  be  promoted  by  allowing  q{ x,  y;t)  ~  qo(t).  When 
(x,  y)  lie  on  an  edge,  then  defining  q(x,  y,  t)  S>  qo(t)  would  pro¬ 
hibit  edge  deterioration.  The  function  q(-)  is  defined  as 


q{x,y\t) 
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SRAD  requires  initialization  of  qo(t),  which  is  determined  by  cal¬ 
culating  the  mean  and  standard  deviation  within  a  homogeneous 
region  where  speckle  is  prevalent. 


3.  THE  SQUEEZE  BOX  FILTER 

The  methods  in  section  2  evaluate  every  sample  in  a  image  and 
adaptively  determines  whether  to  smooth  (locally  average)  or  not. 
The  iterative  filtering  method  we  presented  in  this  section  only 
considers  samples  which  are  outliers  of  some  probability  density 
function  (PDF)  and  applies  local  smoothing  to  these  outliers.  The 
local  extrema  are  considered  outliers  and  are  not  used  in  the  de¬ 
termination  of  the  local  mean.  The  choice  of  the  neighborhood 
A f  is  extremely  important,  since  the  mean  of  some  PDF  is  deter¬ 
mined  by  samples  in  A f.  Each  iteration  produces  a  sequence  with 
locally  reduced  variance.  The  local  extrema  of  the  new  sequence 
are  consider  as  outliers  and  the  process  is  iterated.  The  steps  of  our 
proposed  iterative  method  are  as  follow: 

1.  Each  iteration  i  begins  by  determining  the  set  of  locations 
of  local  maxima  and  local  minima.  The  locations  of  these 
extrema  are  defined  by  the  set 

Me  =  {( n,  m)  |  gi-i(n,  m)  meets  condition  1  or  2  } 

Condition  1:  m)  >  +  k,  m  +  l) 

Condition  2:  g,-i(n,  m)  <  g,_i(n  +  k,  m  +  l) 
where  k,  l  =  —  1  or  1. 

2.  Without  using  the  local  extrema  values,  our  algorithm  re¬ 
places  the  extremum  with  the  local  mean  taken  from  neigh¬ 
boring  samples.  For  all  (n,  m  )  G  Me 

gi{n,m)=-r^ -  ^  fli-i(M)  (25) 


are  shown  along  the  column  at  lateral  position  -15mm  and  at  ax¬ 
ial  distances  40mm,  50mm,  60mm,  70mm,  and  80mm.  Fig.  1(c) 
shows  the  SRAD  result  applied  to  the  simulated  image  in  Fig.  1(b). 
Fig.  1(d)  is  the  resulting  image  of  SBF  applied  to  the  same  image. 
It  can  be  observed  that  the  large  intensity  (white)  disks  on  the  right 
side  of  Fig.  1(d)  is  brighter  and  the  background  class  is  smoother 
with  SBF  restoration. 
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(a)  The  Phantom  (b)  Field  II  simulation 


where  M  is  some  local  neighborhood  of  (n,  m),  |A/”|  is  the 
cardinality  of  set  M,  and  ( n ,  m)  (f:  M . 

3.  If  convergence  is  not  attained,  that  is 

y  m)  —  gi(n,  m)\  >e  (26) 

Vn,m 

for  some  predefined  t  >  0,  then  another  iteration  is  per¬ 
formed.  If  convergence  is  attained,  then  no  further  improve¬ 
ments  can  be  attained  with  this  filtering  method. 

By  removing  outliers  at  each  iteration,  this  method  reduces 
the  local  variance  at  each  pixel.  In  effect,  this  method  produces  a 
convergent  sequence  of  images  by  squeezing  the  stochastically  dis¬ 
tributed  pixel  values  to  a  limiting  value.  Thus,  we  call  this  stochas¬ 
tically  driven  method  the  squeeze  box  filter  (SBF). 

4.  EXPERIMENT  AND  RESULTS 


To  evaluate  the  performance  of  the  proposed  SBF  against  the  other 
methods  in  section  2,  a  phantom,  consisting  of  a  background  class 
with  pixel  values  set  at  one  and  two  other  classes  with  pixel  values 
set  at  either  one  or  ten,  is  used.  Each  nonbackground  class  consist 
of  five  bright  disks  or  five  dark  disks  of  various  diameters  vertically 
aligned. 

The  phantom  and  the  Field  II  simulation  [11],  which  is  the 
image  being  processed  in  this  evaluation,  are  shown  in  Figs.  1(a) 
and  1(b),  resp.  The  simulation  is  constructed  with  the  transducer 
at  the  top  of  the  image.  The  focus  point  of  the  simulation  is  set 
at  70mm  axial  distance  from  the  transducer  and  at  lateral  position 
0mm.  In  this  simulation,  the  spatial  varying  point  spread  functions 


(c)  SRAD  (d)  SBF 

Fig.  1.  (a)  Three  class  phantom,  (b)  Field  II  simulation,  (c)  output 
due  to  SRAD,  and  (d)  output  due  to  SBF. 


To  quantitatively  evaluate  which  despeckling  method  provides 
the  best  improvements  to  an  ultrasound  image,  we  use  a  modified 
Fisher  discriminant  contrast  metric.  This  measure  determines  how 
well  an  algorithm  reduces  variances  in  homogeneous  classes  while 
keeping  the  distinct  classes  well  separated  and  preserving  edges. 


This  metric  is  defined  as 


5.  CONCLUSION 
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where  Ja;s  is  the  resulting  image  due  to  algorithm  alg,  and  \Ck 
denotes  the  number  of  pixels  in  class  Ck  ■  To  avoid  sensitivity  to 
resolution,  we  normalize  the  measure  in  equation  (27)  as 


Q(alg) 


dej  Q(alg) 
Q{id ) 


(30) 


An  overview  of  some  of  the  prominent  speckle  removing  algo¬ 
rithms  for  various  imaging  modalities  is  presented.  These  de- 
speckling  methods  were  used  to  assess  the  performance  of  a  novel 
stochastically  driven  SBF  method.  A  visual  inspection  shows  that 
SRAD  is  better  at  preserving  relevant  point  scatters  than  SBF. 
A  modified  Fisher  discriminant  metric  was  used  to  quantify  the 
contrast  improvement  performance  of  each  despeckling  algorithm. 
The  quantitative  evaluation  using  the  USDSAI  metric  shows  that 
the  SBF  method  performed  the  best  contrast  improvement  to  the 
Field  II  simulated  image  while  reduction  to  intraclass  variance  is 
on  par  with  SRAD. 
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Abstract — A  novel  stochastically  driven  filtering  method  to 
despeckle  B  mode  ultrasound  images  is  presented.  This  method  is 
motivated  by  viewing  the  pixel  values  as  a  stochastic  process  and 
removing  outliers,  where  outliers  are  defined  by  local  extrema. 
These  outliers  are  removed  by  local  averaging.  This  produces 
another  image  with  new  outliers  (local  extrema)  and  the  process 
is  iterated.  With  each  iteration  homogeneous  regions  become 
smoother  while  edges  that  defined  these  regions  are  preserved. 
By  allowing  a  dynamically  varying  window  to  determine  the  local 
mean,  we  achieve  equivalent  results  with  fewer  iterations. 

I.  Introduction 

Speckle  is  a  common  phenomena  found  in  many  imaging 
modalities  such  as  optical  laser,  synthetic  aperture  radar,  and 
ultrasound.  Goodman’s  intuitive  explanation  of  the  appearance 
of  speckle  is  due  to  the  signal  being  composed  of  a  sum 
of  independent  complex  components  caused  by  the  random 
roughness  of  the  object  being  imaged  [1],  [2].  Speckle  can  be 
modeled  as  a  complex  random  walk  where  individual  complex 
components  are  independent  of  each  other  and  the  phase  and 
amplitude  of  each  component  are  also  independent.  Since 
speckle  is  a  common  and  prevalent  phenomena  in  ultrasound 
imaging,  Goodman’s  intuition  and  analysis  is  relevant. 

In  Fig.  1,  we  show  a  block  diagram  of  the  ultrasound  data 
flow  from  the  real-valued  digitized  RF  signal  detected  at  the 
transducer  to  the  intermediate  complex-valued  IQ  data  to  the 
real-valued  B  mode  data.  The  processing  that  occurs  between 
the  RF  and  IQ  data  is  proprietary  to  each  ultrasound  manufac¬ 
ture  and  is  generally  privileged  information.  The  probability 
distribution  of  the  norm  of  the  IQ  to  produce  B  mode  data 
is  considered  in  Goodman’s  analysis  found  in  [2].  Exactly 
applying  Goodman’s  analysis  to  the  IQ  data,  it  can  be  derived 
that  the  samples  of  homogeneous  regions  of  the  B  mode  data 
are  Rayleigh  or  more  generally  Rician  distributed.  As  with  any 
probability  distribution  function  (PDF),  outliers  will  occur  and 
values  close  to  the  mean  occur  more  frequently.  Determining 
which  samples  are  outliers  and  consistently  replacing  these 
outlier  with  a  meaningful  value  is  difficult. 

Numerous  despeckling  methods  indiscriminately  scrutinize 
every  pixel  values  and  rely  on  local  statistics  like  the  local 
mean,  local  variance,  and/or  the  gradient  to  determine  how 
aggressively  a  pixel  should  be  smoothed.  These  local  statistics 
are  usually  determined  by  a  fixed  window  size  usually  of  odd 


dimensions  3  x  3,  5  x  5,  and  so  on.  One  of  the  first  of  these 
filters  is  the  Nagao  and  Matsuyama  filter  proposed  in  [3], 
which  recursively  replaces  each  pixel  value  with  the  mean 
of  the  minimum  variance  segment  from  all  eight  direction 
originating  from  the  pixel  under  scrutiny. 

A  more  profound  and  adaptive  despeckling  method  was 
proposed  by  the  Lee  filter  in  [4],  This  algorithm  adaptively 
chooses  a  weighting  factor  k  between  zero  and  one  so  that 
the  new  pixel  value  J (n,  to)  is  set  at 

J(n,m)  =  p  +  k(J{n,m)  —  /i)  (1) 

where  /i  is  the  mean  in  some  window.  The  weighting  factor 
k  is  adaptively  determined  as 


—  2  K  ' 
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where  <r2  is  the  local  variance  in  the  some  window.  The  noise 
variance  is  denoted  as  er2  and  must  be  known  a  priori.  When 
o2  cr2  ^  0,  the  gain  parameter  is  approximately  one,  in 
which  case  the  filter  in  equation  (1)  is  the  identity  filter.  If  the 
local  variance  o2  is  greater  than  but  nearly  equal  to  the  global 
noise  variance  a2,  then  the  filter  specified  by  equation  (1)  is 
a  local  averaging  low  pass  filter. 

The  Frost  filter  given  in  [5]  consists  of  determining  a  filter 
so  that  in  a  local  homogeneous  region  and  in  the 
presence  of  white  noise,  the  expection  criterion  is  minimized. 
The  term  I (n,  to)  is  a  windowed  weighted  sum  about  some 
constrained  ideal  image  I(n,m .) 

I(n,m)  =  f(n',m')I(n  +  n',m  +  m')  (3) 

n'  ,m'  (zW 

where  W  is  some  odd  dimensional  window. 

The  linear  minimum  mean  squared  error  (LLMSE)  estimate 
of  the  ideal  noise  free  image  given  by  Kuan  et  al.  in  [6]  is 
defined  as 

~  U)2  —  <72 

I(n,m)  =  /-H - ^(J(n,ra)  -  /z)  (4) 

where  p  is  the  local  variance  in  some  window  about  J(n,  to) 
and  lo2  is  the  local  variance  in  some  weighted  window  about 
J(n, to).  The  local  weighted  variance  lo1  is  defined  as 

w2  =  WM1  ^  w{n\m'){J{n  +  n',m  +  rn')-pf  { 5) 

n'  ,m' 


where  y~^  w(n',  m!)  =  1  and  W  is  some  odd  dimensional 

n'  ,m' 

square  window. 

The  adaptive  weighted  median  filter  (AWMF)  in  [7]  models 
the  ultrasound  image  as 

J(n,m)  =  I(n,m)  +  \/l(n,m)r](n,in).  (6) 


From  the  image  model  given  in  equation  (6)  and  when  the 
ideal  image  is  equal  to  a  constant  in  some  local  neighborhood, 
then 

a  =  cav  (7) 

where  a2  and  a2  are  the  local  variance  of  the  observed  image 
and  the  variance  of  the  noise,  resp.  The  AWMF  is  defined  as 

J(n,  to)  =  median^  J(n' ,  to/), . . . ,  J{n\  in')  |  (8) 

w(n'  ,m') 


for  n',m'  in  some  window  about  ( n,m .)  and 

c(T2v/tr'2  +  to'2  1 
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where  round+{-}  means  round  to  the  nearest  non-negative 
integer,  c  is  some  constant,  p  is  the  local  mean,  and  cr2  is  the 
local  variance. 

The  filtering  method  of  Yu  and  Acton  in  [8],  the  speckle 
reducing  anisotropic  diffusion  (SRAD)  method  incorporates 
gradient  information  as  well  as  local  statistics  to  determine 
a  smoothing  kernel.  The  SRAD  algorithm  iteratively  filters  a 
nonzero  valued  image  J(x,y,  0)  =  J(x,y )  according  to 


w(n',m')  =  round+  <  w(0,0) 


and 


9J{x,y;  t ) 
dt 


=  div[c{q)S7  J(x,y;t)] 


dJ (x,  y\  t) 


dn 


=  0 


B(fi) 
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where  n  is  the  outward  normal  vector  to  the  border  B(Q).  The 
diffusion  coefficient  c(-)  is  defined  either  as  the  quotient 


c{q{x,y,t)) 
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9o(f)  +q2(x,y\  t) 
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or  as  the  exponential  function 


c(q(x,  y,  t))  =  exp 


( qo(t )  -q2(x,  y>t)\ 

\  q40(t)  +  q2o(t)  )' 


(13) 


In  both  equations  (12)  and  (13),  if  q(x,y\t )  «  qo(t ),  then 
c(q(x,y,t))  ~  1  and  equation  (10)  is  a  local  smoothing 
operator.  If  q(x,y,t)  qo(t),  then  the  diffusion  coefficient 
is  very  small  and  smoothing  in  a  local  region  around  ( x ,  y) 
is  suppressed.  When  at  time  t,  if  (x,  y)  resides  in  a  homoge¬ 
neous  region,  then  smoothing  can  be  promoted  by  allowing 
q(x,y,t)  ~  qo(t).  When  ( x,y )  lie  on  an  edge,  then  defining 
q{x,y\t)  qo(t)  would  prohibit  edge  deterioration.  The 
function  </(•)  is  defined  as 


q(x,y,  t) 


N 


1  ( wh)2  _  j_  (v2j\2 

2  \  J  J  16  \  J  ) 

(i+a)(^))2 


(14) 


SRAD  requires  initialization  of  qo(t),  which  is  determined 
by  calculating  the  mean  and  standard  deviation  within  a 
homogeneous  region  where  speckle  is  prevalent. 

Determining  the  appropriate  window  to  determine  local 
statistics  is  important  to  these  and  other  algorithms.  We  offer 
a  novel  despeckling  method  that  iteratively  removes  outliers 
by  determining  the  local  mean  and  standard  deviation  from  an 
adaptively  varying  window.  The  adaptively  determined  mean 
is  used  to  replace  the  outlying  values  of  an  B  mode  ultra¬ 
sound  image  causing  homogeneous  regions  to  be  aggressively 
smoothed  while  preservation  of  edges  is  profoundly  respected. 

II.  The  Squeeze  Box  Filter 

We  propose  an  iterative  filtering  method  that  considers  local 
extrema  of  the  B  mode  image  J(n,m .)  as  outliers  and  only 
applies  local  smoothing  to  these  outliers.  The  local  extrema 
are  considered  outliers  and  are  not  used  in  the  determination 
of  the  local  mean.  The  choice  of  the  neighborhood  or  window 
A f  is  extremely  important,  since  the  mean  of  some  PDF  is 
determined  by  samples  in  A f.  An  explanation  of  a  robust 
method  to  determine  this  window  will  be  given  in  section  III. 
Each  iteration  of  our  proposed  method  produces  a  sequence 
with  locally  reduced  variance.  The  local  extrema  of  the  new 
sequence  are  consider  as  outliers  and  the  process  is  iterated. 
The  steps  of  our  proposed  iterative  method  are  as  follow: 

1)  Each  iteration  i  begins  by  determining  the  set  of  loca¬ 
tions  of  local  maxima  and  local  minima.  The  locations 
of  these  extrema  are  defined  by  the  set 

Me  =  {(n,  to)  |  J,;_i  (n.  rri)  meets  condition  1  or  2  } 

Condition  1:  to)  >  Jj_i(n  +  k,  m  +  l) 

Condition  2:  to)  <  Jj_i(n  +  k,m  +  l) 

where  k,  l  =  —1  or  1. 

2)  Without  using  the  local  extrema  values,  our  algorithm 
replaces  the  extremum  with  the  local  mean  taken  from 
neighboring  samples.  For  all  (n,  to)  £  Me 

Ji(n,m)  =  jij-  ^  Ji-i{k,l )  (15) 

1  V  1  (fc.QeW 

where  M  is  some  local  neighborhood  of  (n,  to),  \M\  is 
the  cardinality  of  set  A/”,  and  (n,  m)  ^  M. 

3)  If  convergence  is  not  attained,  that  is 

E  \  Ji-i(n,m)  -  Ji(n,m)\  >  e  (16) 

Vn,m 

for  some  predefined  e  >  0,  then  another  iteration  is 
performed.  If  convergence  is  attained,  then  no  further 
improvements  can  be  attained  with  this  filtering  method. 

By  removing  outliers  at  each  iteration,  this  method  reduces 
the  local  variance  at  each  pixel.  In  effect,  this  method  produces 
a  convergent  sequence  of  images  by  squeezing  the  stochasti¬ 
cally  distributed  pixel  values  to  a  limiting  value.  We  call  this 
stochastically  driven  method  the  squeeze  box  filter  (SBF). 


B  mode 


Fig.  1.  RF  data  to  IQ  data  to  B  mode  image. 


III.  Window  Selection 

Given  £  Me  is  an  extremum  value  at  some 

iteration  i,  some  considerations  for  a  neighborhood  or  window 
M  about  m)  to  determine  the  local  mean  are 

1)  the  samples  of  M  should  be  from  the  same  homogeneous 

region  as  and 

2)  to  get  an  accurate  estimate  of  the  actual  mean  within  a 
homogeneous  region,  the  cardinality  |A/"|  should  be  as 
large  as  possible. 

To  effectively  and  dynamically  determine  the  window  M,  a 
multiplicative  image  model  of  equation  (17)  is  used 

J(n,  to)  =  I(n,  m)r)(n,  m)  (17) 

where  I(n,m)  is  the  speckle  free  image  and  rj(n,m)  is  a 
white  stationary  random  process.  This  model  is  supported  by 
the  empirical  evidence  and  deductive  arguments  given  in  [9], 
[10],  In  a  constant  homogeneous  region  Mc,  I(n,  m)  =  c  ^  0 
for  all  (n,  ?n)  £  Mc ,  the  mean  to  standard  deviation  ratio 
f?Wc(n,  ?ro)  is  independent  of  c,  that  is, 

Rmc  {ni m)  =  (18) 

<JMc 

where  HMC  and  a_\rc  are  the  mean  and  standard  deviation 
within  Mc,  resp.  Let  R  be  the  ratio  defined  by  equation  (18) 
determined  from  a  known  non-zero  constant  region.  At  each 
outlier,  we  dynamically  choose  a  rectangular  window  of  di¬ 
mension  N'  x  M'  so  that  |i?  —  Rnxm\  is  minimized  where 
Rn  x  m  is  the  ratio  of  the  mean  to  standard  deviation  within  a 
rectangular  NxM  window  centered  at  the  outlier  m). 

The  SBF  despeckling  method  using  this  dynamically  varying 
window  (SBF-DVW)  replaces  each  outlier  values  with  the 
mean  determined  from  this  dynamically  varying  window  less 
the  outlier  value. 

IV.  Experiments  and  Results 

To  evaluate  the  performance  of  our  algorithm  we  created 
a  phantom  and  the  Field  II  simulation  [11]  of  the  phantom, 
which  is  the  image  being  processed  in  this  evaluation.  The 
phantom  and  the  Field  II  simulation  are  shown  in  Figs.  2(a) 
and  2(b),  resp.  The  phantom  consist  of  the  background  class 
with  pixel  values  set  at  one.  The  phantom  also  consist  of 
two  other  classes  with  pixel  values  set  at  either  one  or  ten. 
Each  nonbackground  class  consist  of  five  bright  disks  or 
five  dark  disks  of  various  diameters  vertically  aligned.  The 
Field  II  simulation  of  an  actual  B  mode  ultrasound  image 
is  constructed  with  the  transducer  at  the  top  of  the  image. 
The  focus  point  of  the  simulation  is  set  at  70mm  axial 
distance  from  the  transducer  and  at  lateral  position  0mm.  In 
this  simulation,  the  spatial  varying  point  spread  functions  are 
shown  along  the  column  at  lateral  position  -15mm  and  at  axial 
distances  40mm,  50mm,  60mm,  70mm,  and  80mm.  Fig.  2(c) 


shows  the  SRAD  result  applied  to  the  simulated  image  in 
Fig.  2(b).  Fig.  2(d)  is  the  resulting  image  of  SBF-DVW  applied 
to  the  same  image.  The  SBF-DVW  we  implemented  allowed 
for  a  rectangular  window  with  odd  row  and  column  dimensions 
to  vary  N ,  M  =  3,5,...,  41.  It  can  be  observed  that  the 
large  intensity  (white)  disks  aligned  on  the  center  column  of 
Fig.  2(d)  is  brighter  and  the  background  class  is  smoother  with 
SBF-DVW  method.  The  visual  results  of  SBF  using  a  fixed 
9  x  9  is  almost  indistinguishable  from  the  result  of  SBF-DVW. 
Thus  only  the  result  of  SBF-DVW  is  shown  in  Fig.  2. 

To  quantitatively  evaluate  performance  we  use  a  modified 
Fisher  discriminant  contrast  metric.  This  measure  determines 
how  well  an  algorithm  reduces  variances  in  homogeneous 
classes  while  keeping  the  distinct  classes  well  separated  and 
preserving  edges.  This  metric  is  defined  as 


/RA  =  ]7TT  Jaig{n,m), 

1  fc|  (tt,m)6Ct 


2  def  1 

=  i^i 


'y  (Jalgin,  m)  I^Ck)  i 

(n,m)GCfc 


(21) 


where  3aig  is  the  resulting  image  due  to  algorithm  alg,  and 
Ci-  denotes  the  number  of  pixels  in  class  C),  •  To  avoid  sensi¬ 
tivity  to  resolution,  we  normalize  the  measure  in  equation  (19) 


~ ,  ,  def  Q(alg) 

Q(als)  =  QM 


where  .Tu/  =  J.  This  performance  metric  Q(alg)  is  called 
ultrasound  despeckling  assessment  index  (USDSAI).  If  a  de¬ 
speckling  algorithm  produces  classes  that  are  well  separated, 
then  the  numerator  in  equation  (19)  will  be  large.  If  the 
segmentation  or  restoration  algorithm  produces  small  intra¬ 
class  variances,  then  the  denominator  of  equation  (19)  will 
be  small.  An  algorithm  that  attains  a  large  numerator  and  a 
small  denominator  will  yield  a  large  USDSAI  quantity  Q(alg). 
Large  USDSAI  would  indicate  that  alg  produces  desirable 
restoration  or  enhancement  results. 

An  objective  quantitative  comparison  is  performed  by  vary¬ 
ing  the  parameters  over  a  wide  range  of  reasonable  values 
for  each  algorithm  listed  in  Table  I,  the  number  of  iteration 
used  in  SBF  with  a  fixed  9x9  window,  and  SBF-DVW  so 
that  the  USDSAI  value  Q(alg)  is  maximized.  The  USDSAI 
values  Q(alg)  of  each  filter  are  given  in  Table  I.  The  USDSAI 
in  Table  I  are  the  results  attained  by  an  exhaustive  effort  to 
maximize  this  metric.  The  results  of  the  fixed  window  SBF 


(c)  SRAD  (d)  SBF-DVW 


Fig.  2.  (a)  Three  class  phantom,  (b)  Field  II  simulation,  (c)  output  due  to 

SRAD,  and  (d)  output  due  to  SBF-DVW. 

and  SBF-DVW  exceed  the  results  of  the  other  methods.  The 
result  of  SBF  was  attained  by  using  the  values  within  a  9  x  9 
square  window  minus  the  extremum  value  to  determine  the 
local  mean  at  each  extrema.  The  maximum  attained  USDSAI 
of  the  SBF  is  2.1491  and  took  810  iterations  to  attain.  The 
maximum  USDSAI  of  SBF-DVW  is  slightly  better  at  2.1716, 
but  occurred  with  only  143  iterations.  This  result  is  promising 
in  which  we  maybe  able  to  decrease  the  run  time  of  our 
proposed  SBF  algorithm  while  maintaining  excellent  contrast 
improvements. 

V.  Conclusion 

Based  on  the  intuition  that  speckle  is  a  random  stochastic 
process,  we  develop  a  novel  despeckling  method  SBF,  which 


Q(alg) 

F 

Nagao 

1.0315 

I 

Lee 

1.2007 

L 

Frost 

1.0446 

T 

Kuan 

1.0012 

E 

AWMF 

1.0959 

R 

SRAD 

1.8603 

S 

SBF 

2.1491 

SBF-DVW 

2.1716 

TABLE  I 

USDSAI  FOR  THE  VARIOUS  ALGORITHMS  TESTED. 


is  aimed  at  replacing  outliers  with  the  local  mean.  The  SBF 
and  SBF-DVW  methods  attained  better  qualitative  and  quanti¬ 
tative  results  than  the  other  well  known  published  despeckling 
methods.  By  allowing  a  dynamically  varying  window  in  which 
to  determine  the  local  mean,  SBF-DVW  attain  slightly  better 
results  than  SBF  with  a  fixed  window  with  fewer  iterations. 
This  result  is  promising  in  that  we  can  achieve  a  faster 
algorithm  without  compromising  performance. 

Acknowledgment 

This  work  was  supported  by  NIH  NIBIB  grant  EB001826 
and  US  Army  CDMRP  grant  (W81XWH-04- 1-0240). 

References 

[1]  J.  W.  Goodman,  “Statistical  properties  of  laser  speckle  patterns,”  in  Laser 
Speckle  and  Related  Phenomena ,  J.  C.  Dainty,  Ed.  Berlin:  Springer- 
Verlag,  1984,  pp.  9-75. 

[2]  - ,  “Speckle  phenomena  in  optics:  Theory  and  applications  version 

5.0,”  http://www-ee.stanford.edu/~goodman/,  Aug.  2005. 

[3]  M.  Nagao  and  T.  Matsuyama,  “Edge  preserving  smoothing,”  Computer 
Graphics  and  Image  Processing,  vol.  9,  no.  4,  pp.  394—407,  April  1979. 

[4]  J.  S.  Lee,  “Digital  image  enhancement  and  noise  filtering  by  use  of 
local  statistics,”  IEEE  Trans.  Pattern  Anal.  Machine  Intell.,  vol.  PAMI- 
2,  no.  2,  pp.  165-168,  March  1980. 

[5]  V.  S.  Frost,  J.  A.  Stiles,  K.  S.  Shanmugan,  and  J.  C.  Holtzman,  “A 
model  for  radar  images  and  its  application  to  adaptive  digital  filtering 
of  multiplicative  noise,”  IEEE  Trans.  Pattern  Anal.  Machine  Intell.,  vol. 
PAMI-4,  no.  2,  pp.  157-166,  March  1982. 

[6]  D.  T.  Kuan,  A.  A.  Sawchuk,  T.  C.  Strand,  and  P.  Chavel,  “Adaptive 
restoration  of  images  with  speckle,”  IEEE  Trans.  Acoust.,  Speech,  Signal 
Processing,  vol.  ASSP-35,  no.  3,  pp.  373-383,  March  1987. 

[7]  T.  Loupas,  W.  N.  McDicken,  and  P.  L.  Allan,  “An  adaptive  weighted 
median  filter  for  speckle  suppression  in  medical  ultrasonic  images,” 
IEEE  Trans.  Circuits  Syst.,  vol.  36,  no.  1,  pp.  129-135,  January  1989. 

[8]  Y.  Yu  and  S.  T.  Acton,  “Speckle  reducing  anisotropic  diffusion,”  IEEE 
Trans.  Image  Processing,  vol.  11,  no.  11,  pp.  1260-1270,  November 
2002. 

[9]  R.  F.  Wagner,  S.  W.  Smith,  J.  M.  Sandrik,  and  H.  Lopez,  “Statistics  of 
speckle  in  ultrasound  B-scans,”  IEEE  Trans.  Sonics  Ultrason.,  vol.  30, 
no.  3,  pp.  156-163,  May  1983. 

[10]  R.  F.  Wagner,  M.  F.  Insana,  and  S.  W.  Smith,  “Fundamental  correlation 
lengths  of  coherent  speckle  in  medical  ultrasonic  images,”  IEEE  Trans. 
Ultrason.,  Ferroelect.,  Freq.  Contr.,  vol.  35,  no.  1,  pp.  34^14,  Jan.  1988. 

[11]  J.  A.  Jensen  and  N.  B.  Svendsen,  “Calculation  of  pressure  fields  from 
arbitrarily  shaped  apodized  and  excited  ultrasound  transducers,”  IEEE 
Trans.  Ultrason.,  Ferroelect.,  Freq.  Contr.,  vol.  39,  no.  2,  pp.  262-267, 
Feb.  1992. 


